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Abstract 

In  this  thesis,  I  investigated  the  how  the  life  history  characteristics  of  the  clam  Mya 
arenaria  determine  the  population  response  to  chronic  contaminant  exposure.  To 
predict  the  potential  responses  of  a  broadcast-spawning  life  history  such  as  that  of 
M.  arenaria,  I  surveyed  the  literature  on  a  variety  of  bivalve  species.  By  incorporating 
information  on  growth,  survival,  and  reproduction  into  matrix  population  models  I 
could  evaluate  the  relative  contributions  of  these  factors  to  fitness.  For  broadcast- 
spawners,  long  life  is  an  important  factor  enabling  them  to  gamble  on  rare,  large 
recruitment  events.  Another  important  aspect  of  the  broadcast  spawning  strategy  is 
the  possibility  of  high  variation  in  larval  settlement  from  year  to  year.  I  evaluated 
the  role  that  this  variability  plays  using  a  stochastic  matrix  model,  and  showed  that 
it  tends  to  increase  population  growth  because  of  the  larger  size  of  rarer,  successful 
recruitment  events. 

With  an  understanding  of  how  the  life  history  traits  of  M.  arenaria  might  control 
its  responses  to  change  in  the  environment,  I  analyzed  the  vital  rates  of  clams  at 
clean  and  contaminated  sites.  The  effects  of  contaminants  measured  in  the  lab  do 
not  necessarily  predict  population  condition  in  the  field.  Since  surviving  with  a  long 
life  span  contributes  the  most  to  fitness  in  broadcast-spawning  bivalves,  effects  on 
reproductive  output  and  juvenile  survival,  which  are  strong  in  many  lab  studies,  may 
not  necessarily  play  a  large  role  in  field  populations.  The  life  history  of  this  clam,  with 
natural  variation  in  recruitment  from  year  to  year,  further  reduces  the  population 
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dependence  on  high  reproductive  output  and  larval  survival.  The  combination  of 
little  population-level  relevance  of  the  strongest  contaminant  effects,  and  potential 
contaminant  effects  on  very  important  clam  predators,  could  explain  why  populations 
at  contaminated  sites  are  observed  to  be  growing  the  fastest.  The  interaction  of 
contaminant  exposure  and  normal  ecological  processes  determines  the  overall  impact 
on  the  population. 

Thesis  Supervisor;  Judith  E.  McDowell  Capuzzo 
Title:  Senior  Scientist,  WHOI 

Thesis  Supervisor:  Hal  Caswell 
Title;  Senior  Scientist,  WHOI 
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Chapter  1 


Introduction 


1.1  Life  History 

Life  history  traits  are  parts  of  a  suite  of  co-adapted  factors  that  relate  to  survival 
and  reproduction  of  an  organism  (Stearns,  1992).  Trade-offs  required  between  energy 
invested  into  survival  and  reproduction,  for  example,  prevent  organisms  from  having 
the  optimal  strategy  of  maximal  reproductive  success  at  birth.  Although  there  is  some 
plasticity  in  life  history  traits,  the  type  of  life  history  an  organism  has  impacts  how 
it  can  respond  to  change  in  the  environment,  such  as  anthropogenic  contamination. 
For  example,  some  species  have  life  stages  that  are  more  sensitive  to  disturbances 
than  others.  Among  marine  bivalve  species  there  is  a  broad  range  of  life  history 
characteristics.  Some  are  small  and  short-lived,  and  brood  few  offspring,  while  others 
are  longer-lived,  and  broadcast  spawn  many  offspring.  The  latter  life  history  is  the 
most  common  in  marine  invertebrates  in  general. 

1.1.1  Life  History  of  Soft  Shell  Clams 

The  soft-shell  clam,  M.  arenaria,  is  common  and  commercially  important  in  New 
England.  It  has  a  life  cycle  typical  of  marine  benthic  invertebrates,  with  a  sessile 
adult  phase  and  a  dispersive  larval  phase  (Figure  1.1).  The  reproductive  process 
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DEVELOPMENT 


Figure  1.1:  Life  Cycle  of  My  a  arenaria.  Sessile  adult  and  planktonic  larval  phases, 
and  the  processes  which  affect  each  phase,  are  shown. 
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starts  in  late  winter,  when  gametogenesis  begins  with  female  clams  mobilizing  lipid 
reserves  and  depositing  them  into  the  developing  eggs  (Giese,  1959;  Vassalo,  1973). 
The  lipid  pool  in  the  egg  constitutes  the  energy  reserves  that  will  be  available  to 
a  larva  during  its  early  developmental  stages.  Larval  survival  is  dependent  on  the 
maternal  contributions  of  energy  reserves:  bivalve  larvae  from  larger  eggs  (Kraueter 
et  al,  1982)  or  more  lipid-rich  eggs  (Gallager  et  ai,  1986)  have  a  higher  chance  of 
surviving  to  metamorphosis. 

When  water  temperature  rises  in  the  spring,  sperm  and  oocytes  are  shed  simul¬ 
taneously  into  the  water  column,  where  fertilization  takes  place.  After  about  12 
hours,  trochophore  larvae  form,  which  develop  into  veliger  larvae  about  12  hours 
later.  Female  clams  produce  on  the  order  of  10^  6ggs  per  animal  per  spawning  event 
(Brousseau,  1978),  but  larvae  are  subject  to  intense  mortality  during  the  planktonic 
period.  An  approximate  value  for  larval  survival  is  1%  (Thorson,  1966),  but  calcu¬ 
lations  have  shown  that  for  M.  arenaria  this  value  is  essentially  zero  in  some  years 
(Brousseau  et  ai,  1982).  Mortality  is  principally  due  to  eggs  going  unfertilized,  pre¬ 
dation,  starvation  of  larvae,  and  failure  of  larvae  to  find  a  suitable  settling  spot  (Rum- 
rill,  1990).  Interannual  variation  in  these  factors  can  lead  to  variability  in  recruitment 
from  year  to  year. 

Planktonic  larval  development  continues  for  14-21  days  after  fertilization,  where¬ 
upon  the  velum  is  shed  and  larvae  settle  onto  the  substrate.  During  the  post¬ 
settlement  period,  they  can  continue  to  crawl  and  swim  for  2-5  weeks,  before  set¬ 
tling  permanently  (Newell,  1991).  Newly-settled  clams  are  much  more  susceptible  to 
predators  and  environmental  stress  than  adults.  They  grow  rapidly,  however,  and  es¬ 
cape  these  mortality  factors  with  size.  Growth  slows  when  sexual  maturity  is  reached, 
at  about  40  mm  in  size  or  2  years  of  age,  depending  on  water  temperature  and  food 
supply  (Goshima,  1982).  Adult  clams  can  live  longer  than  12  years  (Newell,  1991). 
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1.2  Contaminant  Effects  on  Bivalves 


Many  intertidal  areas  around  the  world  are  heavily  contaminated  with  lipophilic  or¬ 
ganic  compounds  such  as  polychlorinated  biphenyls  (PCBs)  and  polycyclic  aromatic 
hydrocarbons  (PAHs).  Due  to  their  hydrophobicity,  PCBs  and  PAHs  reside  in  the  sed¬ 
iment  and  porewater  and  linger  in  the  environment  even  when  the  source  is  removed 
(Farrington,  1990).  The  long-term  effects  of  such  chemicals  on  the  organisms  exposed 
to  them  are  currently  an  active  area  of  research.  Bivalves  such  as  Mya  arenaria  are 
frequently  used  as  sentinel  organisms  (organisms  that  can  be  easily  monitored  in  the 
environment  for  exposure,  uptake,  and  effects  of  contaminants)  because  they  accu¬ 
mulate  toxins  readily,  process  them  minimally,  and  their  tissues  reflect  only  local 
conditions  because  they  are  sessile  (NRC,  1991). 

1.2.1  Exposure  and  Uptake 

Over  the  life  cycle  of  the  clam,  routes  of  contaminant  exposure  and  ability  to  ab¬ 
sorb  and  store  contaminants  change.  Adults  are  exposed  to  sediment  and  porewa¬ 
ter  contamination  through  the  integument  and  through  the  gut  lining  from  feeding 
on  contaminated  sediment.  Uptake  of  lipophilic  contaminants  is  controlled  by  sim¬ 
ple  equilibrium  partitioning  between  concentrations  in  the  water  and  those  in  the 
hemolymph  and  tissues,  modulated  by  exposure  duration  and  physiological  state  of 
the  animal  (Farrington,  1990).  The  ultimate  fate  of  hydrophobic  contaminants  is 
storage  in  fatty  tissues  such  as  the  gonads  (Moreno  et  al.,  1992).  Total  tissue  lipid 
content  has  been  shown  to  dictate  the  amount  of  hydrocarbon  that  oysters  are  able 
to  accumulate  (Stegeman  and  Teal,  1973). 

Female  clams  mobilize  their  lipid  stores  into  egg  production,  so  some  of  their 
stored  contaminant  load  is  shed  as  they  spawn,  causing  seasonal  differences  in  tissue 
contaminant  concentrations  (Hummel  et  a/.,  1989;  McDowell  Capuzzo  et  al,  1989). 
This  transfer  may  only  be  possible  for  certain  compounds,  though;  lower  molecular 
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weight  PCBs  are  transferred  into  eggs  in  the  blue  mussel,  Mytilus  edulis,  but  higher 
molecular  weight  ones  remain  in  adult  tissues  (McDowell  Capuzzo  et  ai,  1989).  Most 
of  the  contaminant  burden  in  a  larva  is  that  transferred  to  the  egg  by  the  mother, 
since  only  after  oil  spills  or  other  pollutant  input  events  are  hydrophobic  contaminants 
in  the  water  column.  The  fate  of  contaminants  in  larval  tissues  as  they  develop  and 
use  up  lipid  reserves  is  unknown. 

After  settlement,  juvenile  clams  will  start  to  accumulate  lipophilic  contaminants 
from  the  environment  if  they  are  exposed  to  them.  Experiments  have  shown  that 
if  bivalves  are  placed  in  flowing,  clean  water,  they  are  able  to  depurate  their  tissues 
of  low  molecular  weight  oils  in  days  to  weeks  (Widdows  et  al,  1985;  Stegeman  and 
Teal,  1973;  Moreno  et  ai,  1992).  However,  even  low  level  exposures  sustained  over 
time,  the  most  environmentally  relevant  scenario,  cause  concentrations  in  the  animal’s 
tissues  to  persist  (Moreno  et  ai,  1992).  Heavier  oils  and  PCBs  are  more  persistent 
in  tissues.  For  example,  PCB  concentrations  did  not  change  in  mussels  transplanted 
from  a  contaminated  site  to  a  clean  site  (Hummel  et  al,  1989). 

Uptake  of  contaminants  depends  on  the  mixture  of  compounds  present,  and  its 
physical  and  chemical  bioavailablity  to  clams.  Lower  molecular  weight  compounds 
are  volatile  and  more  water  soluble,  and  are  not  extensively  accumulated  in  animal 
tissues,  while  large  compounds  are  not  assimilated  because  of  steric  hindrances  in 
transport  through  membranes  (Farrington  et  ai,  1986).  For  example,  Stegeman  and 
Teal  (1973)  found  that  petroleum  hydrocarbon  fractions  in  oyster  tissues  had  a  greater 
aromatic  content  than  the  oil  to  which  they  were  exposed.  Gardner  and  Pruell  (1988) 
found  that  PCB  congeners  in  M.  arenaria  tissues  were  of  lighter  molecular  weight 
than  those  found  in  sediments.  Further,  concentrations  in  clam  tissues  were  2-3  times 
higher  those  found  in  sediments,  demonstrating  bioaccumulation. 
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1.2.2  Metabolism  of  Lipophilic  Contaminants 

Excretion  of  toxins  and  potential  for  damage  depend  on  the  extent  of  metabolic  trans¬ 
formation  of  the  compounds.  Metabolism  of  lipophilic  contaminants,  as  is  known  in 
vertebrates,  is  a  two-step  process.  First,  compounds  undergo  oxidation,  reduction, 
or  hydrolysis  reactions  (phase  I)  to  expose  functional  groups  necessary  for  the  con¬ 
jugation  or  synthesis  reactions  which  constitute  phase  II.  Initially,  researchers  had 
difficulty  demonstrating  ability  of  bivalves  to  transform  xenobiotics  but  there  is  now 
evidence  that  they  at  least  possess  the  appropriate  enzyme  systems.  Activity  of  cy¬ 
tochrome  P-450,  a  well-known  phase  I  enzyme  in  vertebrates,  has  been  detected  in  the 
digestive  gland  of  mussels  and  clams  (Stegeman,  1985;  Livingstone  and  Farrar,  1985), 
although  levels  of  activity  are  an  order  of  magnitude  lower  than  in  mammals  (Ribera 
et  al,  1989).  Activity  of  other  phase  I  enzymes  has  also  been  reported,  particularly 
benzo[a]pyrene  hydroxylase  (BPH:  Moore  et  al.,  1980;  Anderson,  1985).  These  re¬ 
sponses  have  all  shown  some  level  of  induceability  in  the  lab,  but  variability  is  high. 
In  the  field,  no  clear  dose  response  to  environmental  exposure  has  been  observed, 
limiting  the  utility  of  measuring  phase  I  enzymes  as  biomarkers  for  contaminant  ex¬ 
posure. 

Most  of  the  cellular  damage  from  PCBs  and  PAHs  results  from  the  highly  re¬ 
active  intermediate  compounds  generated  by  the  phase  I  reactions.  Evidence  of 
these  reactive  intermediate  metabolites  would  demonstrate  the  extent  of  xenobiotic 
metabolism  occurring  in  bivalves,  but  few  researchers  have  studied  it.  Anderson 
(1985)  found  the  highly  mutagenic  7,8-diol  BPH  product  in  hard  clams  {Mercenaria 
mercenaria)  exposed  to  Aroclor  1254,  but  it  was  also  found  in  controls.  Mytilus  edulis 
converted  about  17%  of  an  aromatic  amine  body  burden  to  N-acyl  derivatives,  which 
can  form  mutagenic  DNA  adducts  (Knezovich  et  al,  1988).  Kurelec  and  Britvic 
(1985)  showed  that  digestive  gland  microsomes  of  M.  galloprovincialis  also  activated 
aromatic  amines,  but  did  not  activate  PAHs.  The  reaction  was  non-inducible,  in¬ 
dicating  that  its  normal  physiological  function  may  not  relate  to  xenobiotics.  It  is 
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hypothesized  that  the  lower  toxicity  levels  often  reported  for  bivalves  may  be  due  to 
their  lower  capacity  to  convert  contaminants  to  damaging,  highly  reactive  intermedi¬ 
ate  forms  (Hale,  1989). 

Reactive  intermediates  generated  by  phase  I  reactions  are  scavenged  by  phase 
II  enzymes,  which  catalyze  reactions  yielding  excretable  hydrophilic  molecules.  An 
example  of  a  phase  II  enzyme  is  glutathione-S-transferase  (GSH),  which  scavenges 
damaging  free  radicals  and  forms  hydrophilic  glutathione  conjugate  products.  A 
freshwater  mussel  showed  an  order  of  magnitude  greater  level  of  activity  of  GSH  af¬ 
ter  exposure  to  dieldrin  or  lindane.  (No  metabolites  of  these  organochlorine  toxins 
were  detected,  however,  suggesting  that  they  are  not  actually  being  biotransformed: 
Boryslawsyj  et  al,  1988).  Ribera  et  al.  (1989)  found  that  due  to  higher  levels  of 
free  radical  scavengers  (GSH,  retinol,  and  tocopherol)  and  fewer  activated  molecules 
because  of  lower  P-450  activity,  M.  galloprovincialis  had  the  same  rate  of  lipid  per¬ 
oxidation  as  rats.  Regrettably,  such  measures  of  contaminant  effects,  using  rates  of 
xenobiotic  metabolism,  have  not  yet  characterized  bivalve  responses  to  environmental 
gradients  sufficiently  for  use  in  monitoring  programs.  Levels  of  heat  shock  proteins 
and  DNA  strand  breakage  are  two  other  potential  biochemical  biomarkers  for  con¬ 
taminant  exposure.  Results  from  preliminary  studies  of  these  responses  suffer  from 
incompletely-determined  methodologies  and  poor  controls,  but  may  prove  useful  after 
further  research  (Clayton,  1996). 

1.2.3  Physiological  Effects 

Correlations  between  contaminant  exposure  and  various  pathologies  and  physiological 
impairments  have  been  more  clearly  documented.  For  example,  the  presence  of  a  suite 
of  histopathological  disorders  was  sufficient  to  discriminate  between  sites  grouped  into 
three  levels  of  contamination  (Moore  et  al,  1996).  Gardner  et  al  (1991)  found  that 
laboratory  and  field  exposure  of  oysters  to  sediment  laden  with  PCBs,  PAHs  and 
heavy  metals  (including  known  carcinogens  and  tumor  promoters)  caused  neoplas- 
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tic  disorders  after  one  month,  while  no  lesions  appeared  in  control  animals.  Clams 
in  contaminated  sites  may  be  more  susceptible  to  disease  (Mix,  1988)  or  parasites 
(Gardner  and  Pruell,  1988)  than  animals  in  clean  sites. 

Effects  on  reproduction  due  to  pollution  exposure  are  linked  to  changes  in  en¬ 
ergy  storage  at  the  cellular  and  organismal  level.  Energy  budget  measurements  in 
M.  edulis  showed  decreasing  scope  for  growth  with  increasing  PAH  levels,  both  in 
the  field  and  in  mesocosm  experiments  (Widdows  and  Johnson,  1988).  Oysters  in 
areas  of  high  PAH  had  lower  lipid  levels  (Bender  et  ai,  1988),  and  exposure  of  blue 
mussels  to  hydrocarbons  caused  a  reduction  in  lipid  storage  levels  including  oocyte 
atresia  (breakdown)  and  resorption  (Lowe  and  Pipe,  1987).  Lowe  and  Pipe  (1987) 
also  showed  that  effects  of  contaminants  are  modulated  by  the  energy  reserves  of  the 
exposed  animal.  In  their  study,  pre-spawning  mussels  had  a  mortality  rate  of  30% 
after  high  doses  of  PAH,  while  post-spawn  mussels  had  a  mortality  rate  of  70%.  They 
hypothesized  that  the  pre-spawn  animals  were  able  to  survive  better  by  remobiliz¬ 
ing  energy  stored  in  gametes  to  meet  metabolic  demands,  whereas  this  energy  was 
no  longer  available  to  the  post-spawn  animals.  Oocyte  resorption  is  not  a  conscious 
choice  of  the  organism,  since  it  is  probably  due  to  xenobiotic-induced  destabilization 
of  lysosomal  membranes.  Once  membranes  are  compromised,  digestive  enzymes  con¬ 
tained  in  the  lysosome  attack  the  cell  itself.  Lipids  are  resorbed  along  with  other 
cellular  components  and  can  then  be  redistributed  to  other  tissues  (Moore,  1988). 

Studies  determining  LC50  values  for  toxicity  tests  have  shown  that  environmental 
levels  of  some  heavy  metals  (Martin  et  ai,  1981),  oils  (Bryne  and  Calder,  1977),  and 
tri-butyl-tin-oxide  (Beaumont  and  Budd,  1984)  will  kill  bivalve  larvae,  but  sublethal 
effects  are  very  poorly  studied.  Bryne  and  Calder  (1977)  showed  that  high  concen¬ 
trations  of  oils  caused  M.  mercenaria  embryos  to  cease  dividing  and  fall  apart  after 
three  hours.  This  study  also  showed  that  if  there  were  a  decrease  in  survival  of  larvae 
in  a  treatment,  there  would  be  slower  growth  in  the  survivors.  A  study  on  larval 
viability  in  hard  clams  from  five  sites  in  Connecticut  demonstrated  increased  levels  of 
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embryo  abnormalities  such  as  chromosome  stickiness  at  contaminated  sites  (Stiles  et 
ai,  1991).  The  only  other  study  of  sublethal  effects  on  bivalve  larvae  demonstrated 
that  exposure  to  contaminated  sediments  causes  failure  to  metamorphose  in  larvae  of 
Cmssostrea  gigas,  the  Japanese  oyster  (Phelps  and  Warner,  1990).  Despite  a  lack  of 
data  on  the  mechanisms,  larvae  are  generally  considered  the  most  sensitive  life  cycle 
stage  to  contaminant  exposure. 

Interactions  between  the  physiological  problems  discussed  above  remain  to  be 
determined,  despite  evidence  that  they  are  important.  Kluytmaus  et  al.  (1988) 
showed  that  the  effects  of  cadmium  exposure  on  M.  edulis  were  both  an  increase  in 
pre-spawn  mortality  and  an  increase  in  the  number  of  spawning  individuals.  Although 
it  is  possible  that  these  effects  cancel  one  another  out,  the  relative  magnitudes  of  the 
changes  need  to  be  analyzed  before  any  such  generalizations  are  made.  Further,  some 
potential  contaminant  effects  may  seem  counter-intuitive.  For  example,  reduced  egg 
production  may  actually  lead  to  increased  fertilization  success,  if  sperm  density  is  a 
limiting  factor  (Levitan  and  Peterson,  1995).  This  type  of  information  reinforces  the 
argument  for  studying  contaminant  effects  in  an  integrated  manner. 

1.2.4  Effects  on  Population  Processes 

Contaminant  effects  may  impact  not  only  the  individual,  but  also  the  population. 
Population  growth,  at  the  most  basic  level,  is  defined  as  the  number  of  births  minus 
the  deaths.  Reduction  in  the  number  of  births  is  a  proven  consequence  of  contaminant 
exposure,  as  discussed  above.  The  impact  on  the  population  will  depend  on  the  extent 
of  reduction  and  the  life  history  strategy  of  the  organism.  For  example,  the  broadcast 
spawning  bivalve  life  history  strategy  with  high  reproductive  output  and  a  planktonic 
larval  phase  may  not  be  sensitive  to  reductions  in  births,  since  so  few  larvae  survive 
to  metamorphosis. 

Deaths  can  be  increased  by  contaminant  exposure,  although  direct  mortality  from 
contaminant  exposure  is  rare,  except  in  cases  of  sudden  high  exposure,  such  as  an  oil 
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spill  at  a  previously  uninsulted  site.  Barry  and  Yevich  (1975)  and  Appledoorn  (1983) 
report  cases  of  immediate  and  massive  mortality  of  M.  arenaria  after  oil  spills.  Dow 
(1975)  also  found  that  if  clams  were  transplanted  from  a  clean  to  a  contaminated  site, 
mortality  was  60%  higher  than  control  rates.  Timing  of  the  exposure  to  chemical 
stress,  with  respect  to  reproductive  condition  and  energy  reserve  levels,  can  modulate 
ability  to  survive  the  insult  (Lowe  and  Pipe,  1987).  Increased  susceptibility  to  disease 
and  environmental  stresses,  such  as  unusually  cold  winters,  are  more  likely  to  cause 
increased  mortality  rates  where  chronic  contamination  is  present. 

These  studies  do  not  address  the  more  common  situation  where  bivalves  settle 
and  grow  to  maturity  in  a  chronically  contaminated  environment.  In  this  case,  either 
genetic  selection  for  contaminant-resistance  at  the  metamorphosis  stage,  or  physiolog¬ 
ical  adaptation  to  the  habitat  may  play  a  role  in  the  health  of  the  population.  Both 
of  these  processes  are  possible.  For  example,  post-settlement  selection  for  salinity 
tolerance  has  been  shown  to  take  place  in  populations  of  Mytilus  edulis  (Hilbish  and 
Koehn,  1985a  and  b).  Regardless  of  the  process  generating  the  dominant  phenotypes, 
differences  in  individual  physiology  between  clean  and  contaminated  sites  can  lead  to 
different  population  characteristics. 

Slower  growth  rates  are  reported  commonly  in  hydrocarbon  effects  studies  on 
M.  arenaria  (Dow,  1975;  Appledoorn,  1981).  The  impacts  of  slower  growth  and 
smaller  size-at-age  relationships  on  population  reproductive  output  over  time  are 
several.  The  onset  of  sexual  maturity  is  more  size-related  than  age-related  in  bivalves 
{e.g.  Nakaoka,  1993;  Newell,  1991),  so  age  at  first  reproduction  would  increase  with 
slower  growth.  Maturation  time  was  found  to  be  one  of  the  more  important  factors 
in  population  growth  rate  differences  between  treatments  in  a  life  table  response 
experiment  on  polychetes  in  hydrocarbon,  sewage,  and  algal  enrichments  (Levin  et 
al,  1996).  Whether  this  effect  is  as  important  in  bivalves,  which  have  much  higher 
lifetime  reproductive  output  and  more  variable  larval  survival,  remains  to  be  explored 
with  modeling  techniques. 
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Since  reproductive  output  increases  with  adult  size  in  most  bivalves  {e.g.,  Peter¬ 
son,  1986),  slower  growth  rates  would  result  in  a  reduced  individual  reproductive  out¬ 
put.  Peterson  also  showed,  as  have  many  others,  that  gonad  volume  increases  faster 
than  shell  size,  so  the  affect  on  reproductive  output  would  be  larger  than  expected 
from  shell  length  differences.  Sensitivities  of  population  growth  rate  to  reproduction 
in  M.  arenaria  depend  on  whether  the  population  is  growing  or  declining.  The  rel¬ 
ative  importance  of  reproduction  in  the  larger-sized  clams  increases  with  decreases 
in  population  growth  rate  (Brousseau  and  Baglivo,  1984).  This  is  due  to  the  higher 
reproductive  contributions  of  larger  clams. 

Slower  growth  rates  can  also  lead  to  increased  predation  risk,  as  many  bivalves 
attain  a  size  refuge  from  some  predators.  For  example,  Nakaoka  (1993)  showed  that 
only  small  Yoldia  notabilis  (an  infaunal  clam)  were  susceptible  to  predation  from  a 
crab.  Blundon  and  Kennedy  (1982)  demonstrated  that  larger  M.  arenaria  obtain  a 
refuge  from  blue  crab  predation  because  they  can  burrow  deeper  than  smaller  clams. 
Spending  longer  at  smaller  sizes  would  expose  clams  to  predation  risk  for  longer, 
reducing  survival  to  larger,  more  fecund,  size. 

The  interactions  between  different  effects  at  different  stages  of  the  life  cycle  and 
their  relative  importance  are  only  starting  to  be  unravelled.  A  few  studies  have 
shown  mixed  responses  to  contamination,  such  as  reduced  egg  numbers  but  increased 
spawning  frequency,  after  cadmuim  exposure  in  Mytilus  edulis  (Kluytmans  et  ai, 
1988).  The  impact  on  the  population  would  depend  on  the  roles  of  survival  and 
reproduction  in  the  life  history  strategy  of  the  organism.  Investigation  of  population 
level  contaminant  effects  has  barely  begun. 

1.3  Matrix  Population  Models 

Mathematical  models  have  proven  to  be  powerful  tools  in  investigating  complex  bio¬ 
logical  processes,  either  as  purely  theoretical  tools  or  as  methods  for  analyzing  exper- 
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imental  data.  For  example,  Levin  et  al.  (1994)  investigated  population  consequences 
of  sewage  enrichment  on  a  polychaete.  A  demographic  model  of  M.  arenaria  has 
clearly  linked  the  disease  state  of  individuals  to  effects  at  the  population  level  (Wein¬ 
berg  et  al.,  1997).  A  series  of  structured  models  on  contaminant  effects  on  populations 
of  Daphnia  have  shown  how  the  lipid  pools  of  individuals  mediate  their  response  to 
lipophilic  contaminants,  and  how  the  structure  of  populations  made  up  of  individuals 
in  differing  physiological  states  are  altered  by  contaminant  exposure  (Lassiter  and 
Hallam,  1990;  Hallam  et  al,  1990;  Hallam  et  al,  1993).  In  both  types  of  studies, 
analysis  made  possible  by  modelling  techniques  has  given  results  that  extend  beyond 
isolated  experimental  results. 

The  impact  of  contaminants  must  be  evaluated  at  each  phase  of  the  life  cycle  to 
understand  how  sublethal  effects  on  individuals  can  be  translated  to  the  population 
level.  Stage-classified  matrix  population  models  integrate  information  on  all  life  cycle 
stages  particularly  clearly  (Lefkovich,  1965).  Analysis  of  the  model  provides  informa¬ 
tion  on  the  growth  rate  of  the  population,  the  impact  on  it  of  certain  perturbations  in 
vital  rates,  and  how  the  vital  rates  interrelate  (Caswell,  1989).  The  model  also  allows 
assessment  of  the  biological  significance  of  statistically  significant  changes  (Weinberg 
et  al,  1986).  While  other  types  of  mathematical  models  can  also  do  this,  the  amount 
of  data  they  require  is  a  deterrent  from  their  application.  Detailed,  species-specific 
structured  models  have  elucidated  population  processes  for  Mytilus  edulis  (Ross  and 
Nesbit,  1990)  and  cladocerans  (Nisbet  et  al,  1989;  Hallam  et  al,  1990),  but  they 
require  many  physiological  and  environmental  parameters.  Matrix  models  are  simple 
to  analyze,  and  they  only  require  data  on  growth,  survival,  and  fecundity,  while  still 
providing  insight  into  fundamental  population  processes.  This  type  of  model  has  al¬ 
ready  been  used  by  several  researchers  to  study  bivalves  (Nakaoka,  1993;  Weinberg  et 
al,  1989;  Weinberg  et  al,  1996;  Brousseau  et  al,  1982;  Brousseau  and  Baglivo,  1984, 
1988;  Malinowski  and  Whitlatch,  1988). 

A  special  application  of  the  matrix  model,  called  life  table  response  experiments 
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(LTREs),  has  been  used  to  analyze  effects  of  contaminants  on  populations  of  some  an¬ 
imals  [e.g.  references  in  Caswell,  1989  and  Caswell,  1996).  LTRE  analysis  interprets 
changes  in  the  the  vital  rates  as  the  response  variable  to  changes  in  environmental 
factors.  For  the  polychaete  Streblospio  benedicti,  hydrocarbon  exposure  was  found 
to  reduce  fertility  and  increase  age  at  first  reproduction  slightly,  leading  to  a  small 
change  in  population  growth  rate.  Sewage  or  algal  enrichment  led  to  decreased  sur¬ 
vival,  but  reduced  age  at  first  reproduction  by  a  factor  of  three  or  two,  the  net  effect 
being  a  large  increase  in  population  growth  rate  (Levin  et  ai,  1994).  Other  LTRE 
studies  have  shown  that  the  most  obvious  effect  of  a  toxicant  on  the  vital  rates  is 
seldom  the  source  of  that  toxicant’s  effect  on  the  population  (Caswell,  1996). 

1.4  Research  Question 

I  investigated  the  population  response  of  M.  arenaria  to  chronic  contaminant  expo¬ 
sure.  I  approached  the  problem  by  considering  the  life  history  of  this  clam  and  how 
the  interaction  of  its  life  history  traits  with  environmental  changes  would  change  pop¬ 
ulation  characteristics.  To  be  able  to  predict  the  potential  responses  of  a  broadcast¬ 
spawning  life  history,  I  surveyed  the  literature  on  a  variety  of  bivalve  species  and 
incorporated  information  on  growth,  survival,  and  reproduction  into  matrix  popu¬ 
lation  models.  This  information  allowed  evaluation  of  the  relative  contributions  of 
these  factors  to  fitness  (Chapter  2).  The  broadcast  spawning  strategy  results  in  the 
possibility  of  high  varation  in  larval  settlement  from  year  to  year.  I  evaluated  the 
role  that  this  variabiliy  plays  using  a  stochastic  matrix  model  (Chapter  3).  With 
an  understanding  of  how  the  life  history  traits  of  M.  arenaria  might  control  its  re- 
sposes  to  change  in  the  environment,  I  analyzed  the  vital  rates  of  clams  at  clean  and 
contaminated  sites  to  test  the  hypothesis  that  population  growth  rates  at  sites  with 
contaminants  would  be  lower  than  those  at  uncontaminated  sites  (Chapter  4). 
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1.5  Life  cycle-environment  interactions  determine 


observed  contaminant  effects 

Results  from  this  thesis  show  that  the  effects  of  contaminants  measured  in  the  lab  do 
not  necessarily  predict  population  condition  in  the  field.  Since  surviving  with  a  long 
life  span  contributes  the  most  to  fitness  in  broadcast-spawning  bivalves  (Chapter  2), 
effects  on  reproductive  output  and  juvenile  survival,  which  are  strong  in  many  lab 
studies,  will  not  necessarily  play  a  large  role  in  field  populations.  The  life  history  of 
this  clam,  with  natural  variation  in  recruitment  from  year  to  year,  further  reduces 
the  population  dependence  on  high  reproductive  output  and  larval  survival  (Chapter 
3).  There  was  no  evidence  that  physiological  changes  in  clams  due  to  contamination 
led  to  the  observed  population  processes  for  M.  arenaria  at  our  study  sites.  I  propose 
that  this  is  due  to  the  combination  of  little  population-level  relevence  of  the  strongest 
known  contaminant  effects,  and  potential  contaminant  effects  on  very  important  clam 
predators.  The  interaction  of  contaminant  exposure  and  normal  ecological  processes 
determines  the  overall  impact  on  the  population  (Chapter  4). 
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Chapter  2 


Contributions  of  growth,  survival, 
and  reproduction  to  fitness  in 
brooding  and  broadcast  spawning 
marine  bivalves 

2.1  Introduction 

Marine  bivalves  are  a  highly  successful  group  of  invertebrates,  found  from  the  high 
intertidal  to  the  deep  sea,  throughout  the  world  ocean.  Their  life  styles  in  these 
varied  environments  include  differences  in  life  history  traits,  such  as  age  at  maturity, 
reproductive  life  span,  reproductive  output,  and  growth  rates.  However,  a  common 
phylogenetic  history  limits  the  extent  of  possible  adaptations  observed  in  bivalves. 
For  example,  the  amount  of  energy  invested  in  eggs  is  limited  by  the  amount  of  space 
available  for  egg  storage  inside  the  shell.  The  diversity  in  life  history  strategies  exhib¬ 
ited  by  this  single  class  makes  it  possible  to  examine  the  demographic  consequences  of 
different  sets  of  life  history  traits.  Anlyzing  only  bivlaves  controls  for  differences  that 
might  be  due  to  phylogeny,  which  would  confound  a  study  comparing,  for  example. 
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bivalves  and  asteroids. 

While  the  comparative  study  of  life  history  strategies  has  long  been  recognized  as 
worthwhile  (Cole,  1954),  it  is  data  intensive  and  conclusive  results  have  been  diffi¬ 
cult  to  obtain.  Studies  on  groups  of  mammals  (Promislow  and  Harvey,  1990),  birds 
(Murphy,  1989),  and  freshwater  bivalves  (Way,  1988)  have  compared  isolated  life  his¬ 
tory  traits,  such  as  mortality  rates  or  energy  invested  in  eggs,  by  standard  statistical 
methods.  The  majority  of  work  incorporating  multiple  life  history  elements  on  more 
than  two  species  at  once  has  been  done  on  plants  [e.g.  Grime,  1977;  Silvertown  et 
al,  1993),  with  the  exception  of  one  study  on  echinoderms  (Ebert,  1996). 

One  method  of  comparing  life  histories  is  to  group  species  according  to  the  envi¬ 
ronment  in  which  they  live,  which  assumes  that  similar  environments  exert  similar 
selective  forces  (Southwood,  1977).  This  method  has  been  applied  to  many  species 
of  plants,  which  are  categorized  as  competitors  (C),  stress-tolerators  (S),  or  ruderals 
(weedy  species)  (R)  by  the  levels  of  particular  traits  that  they  display  (Grime,  1977). 
Amounts  of  stress  and  disturbance  in  an  environment  are  hypothesized  to  determine 
where  each  of  the  three  types  of  life  histories  are  found.  A  triangular  plot  (Figure  2.1) 
with  C,  S,  and  R  axes  emphasizes  the  trade-offs  between  these  strategies.  Groups 
of  plants  from  similar  environments  tend  to  fall  in  particular  regions  on  this  plot, 
demonstrating  that  they  share  many  life  history  traits. 

This  descriptive  method  and  quantitative  demographic  methods  have  recently 
been  reconciled  using  elasticity  analysis  of  matrix  population  models  (Silvertown  et 
al,  1993;  but  see  Shea  et  al,  1994,  for  caveats).  In  a  matrix  population  model,  the 
life  cycle  of  the  organism  is  divided  into  stages  (Figure  2.2),  and  transitions  between 
these  stages  due  to  growth  (G),  survival  (P),  and  reproduction  (F)  are  estimated. 
These  rates  of  transition  are  the  projection  matrix  elements.  The  dominant  eigen¬ 
value,  A,  of  the  projection  matrix,  A,  is  the  asymptotic  population  growth  rate  and 
measures  average  fitness  of  individuals  in  the  population  (Caswell,  1989).  Sensitivity 
of  A  to  change  in  each  matrix  element  can  be  calculated.  Elasticities,  which  are  the 
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Figure  2.1;  Grime’s  Triangle.  The  axes  of  the  triangle  are  the  levels  of  competetive 
strength  (C),  stress-tolerance  (S),  and  ruderality  or  weediness  (R)  found  in  plant 
life  histories.  Circles  show  the  region  of  GSR  parameter  space  occupied  by  various 
groups  of  plants:  a)  annual  herbs;  b)  biennial  herbs;  c)  lichens;  and  d)  trees  and 
shrubs.  Redrawn  from  Grime  et  a/.,  1977. 
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proportional  sensitivities  of  A  to  changes  in  individual  matrix  elements,  can  also  be 
calculated,  and  they  always  sum  to  one  (de  Kroon  et  al,  1986).  Elasticity  analy¬ 
sis  thus  measures  the  relative  contributions  of  different  matrix  elements  to  fitness 
(Caswell,  1989).  By  summing  the  elasticities  of  each  of  the  three  types  of  elements, 
one  can  identify  the  relative  importance  of  growth,  survival,  and  reproduction  to  fit¬ 
ness  of  individuals  in  a  population.  By  equating  the  “competitors”  of  Grime’s  scheme 
to  growth-maximizers,  the  “stress-tolerators”  to  survival-maximizers,  and  “ruderals” 
or  fugitive  species  to  reproduction-maximizers,  the  descriptive  CSR  scheme  can  be 
interpreted  as  a  quantitative  GPF  scheme  (Silvertown  et  al,  1992;  Silvertown  et  al, 
1993).  Plotting  the  sums  of  each  group  of  elasticities  in  the  same  type  of  triangular 
diagram  (Figure  2.3)  used  by  Grime  shows  the  general  pattern  of  importance  of  the 
three  categories  to  fitness  in  each  species. 

A  problem  with  summing  elasticities  to  assess  contributions  to  fitness  is  that 
although  Pi  are  called  “survival”  terms,  survival  is  included  in  the  calculation  of  Pi, 
Fi,  and  Gi  terms.  Differences  in  the  importance  of  survival  for  different  species  are 
confounded  by  simply  summing  types  of  paramerers  (Enright  et  al,  1995).  This  is 
one  of  the  reasons  that  Silvertown ’s  method  of  comparing  elasticities  is  unsatisfying 
to  some  authors  (Shea  et  al.,  1994).  When  matrix  elements  are  calculated  from  terms 
specifying  survival  separately  from  growth,  as  we  have  done  here,  the  elasticities  to 
these  lower-level  parameters  can  be  calculated  (Caswell,  1989).  This  technique  is  used 
rarely  in  matrix  modelling  {e.g.  Brault  and  Caswell,  1993;  Levin  et  al.,  1994),  and 
has  not  been  applied  in  comparative  demography  because  the  elasticities  to  lower- 
level  parameters  do  not  necessarily  sum  to  one  (Caswell,  1989).  Here,  we  calculate 
elasticities  to  lower-level  terms  rrii,  Oi,  and  7ij  used  in  calculating  Gi,  Pi,  and  Fi  and 
investigate  the  utility  of  the  method  in  comparative  demographic  analysis. 

The  broader  goal  of  this  paper  was  to  determine  if  different  types  of  bivalve  life 
histories  have  fundamentally  different  patterns  of  elasticities,  and  thus,  different  con¬ 
tributions  to  fitness  by  the  processes  represented  by  matrix  elements.  Elasticity  anal- 
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P1  P2  P3 


Figure  2.2;  Life  cycle  diagram.  Circles  represent  stages  and  arrows  represent  transi¬ 
tions  between  them  over  one  time  step.  Parameters  are  elements  of  projection  matrix, 

A. 
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Figure  2.3:  Silvertown’s  Triangle.  Axes  are  the  sums  of  elasticities  to  G  (growth), 
P  (stasis),  and  F  (reproduction)  parameters.  The  circles  show  regions  of  parameter 
space  occupied  by  a)  semelparous  herbs;  b)iteroparous  herbs  and  c)  woody  plants. 
Redrawn  from  Silvertown,  1992. 
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ysis  was  performed  on  simple  stage-classified  matrix  population  models  for  ten  bivalve 
species,  over  a  broad  range  of  lifestyles,  parameterized  with  data  from  the  literature. 
The  discussion  of  our  results  focuses  on  the  dichotomy  between  reproductive  modes 
with  broadcast  spawning  and  those  with  brooded  larvae.  Thorson  (1950)  first  grouped 
marine  invertebrates  into  these  categories.  The  former  group  releases  gametes  into 
the  water  column  where  they  develop,  while  the  latter  retains  the  embryos  within  the 
shell  for  part  or  all  of  the  development  time,  releasing  benthic  larvae  or  crawl-away 
juveniles.  Many  authors  have  noted  that  brooders  and  broadcasters  differ  most  in 
juvenile  survival  (higher  for  brooders),  fecundity  (lower  for  brooders),  and  adult  body 
size  (smaller  for  brooders),  but  differences  in  other  factors,  such  as  life  span,  have  not 
been  as  well  studied. 

2.2  Methods 

2.2.1  The  Model 

We  are  interested  in  the  relative  importances  of  the  demographic  processes  of  growth, 
survival,  and  reproduction,  to  fitness  in  marine  bivalves.  We  approached  this  problem 
by  creating  a  simple  demographic  model  in  which  these  processes  could  be  included, 
and  from  which  fitness  could  be  calculated.  The  model  had  to  be  complex  enough 
to  include  the  differences  between  the  species,  simple  enough  that  it  was  possible  to 
parameterize  from  available  data,  and  have  a  form  by  which  all  the  life  histories  could 
be  described.  In  most  cases,  a  bivalve’s  life  span  can  be  divided  up  by  events  such 
as  the  onset  of  sexual  maturity,  and  differences  in  reproductive  output  and  mortality 
with  size.  We  chose  to  use  a  stage-classified  model,  which  divided  the  life  cycle  of 
the  bivalve  into  three  periods  of  time  corresponding  to  periods  with  different  rates  of 
survival  and  reproduction.  The  length  of  time  spent  in  each  stage  depended  on  the 
life  span  of  the  organism  and  the  timing  of  the  life  cycle  events.  More  specifically,  we 
defined  stage  one  as  beginning  at  one  year  from  fertilization.  Stage  two  began  when 
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the  first  of  the  major  life  history  events  occurred,  such  as  a  substantial  change  in 
survival.  Stage  three  began  when  another  of  these  events  occurred,  such  as  attaining 
sexual  maturity.  Bivalves  survived  and  remained  in  stage  three  with  a  probability 
that  was  dependent  on  their  life  span. 

The  equations  describing  the  contribution  of  each  stage  to  the  others  are: 

n(t  +  l)  =  An(t)  (2.1) 


where 

Pi  F2  P3 

A  =  Gi  P2  0 

0  G2  P3 

and  n  is  a  vector  of  stage  abundances.  The  projection  interval  used  was  one  year. 
A  set  of  possible  transitions  of  this  model  are  illustrated  in  a  generalized  life  cycle 
diagram  (Figure  2.2). 

Parameters  were  estimated  from  the  variety  of  types  of  data  that  was  found, 
ranging  from  projection  matrices  to  size-frequency  histograms,  which  are  described 
in  detail  (below)  for  each  species.  Methods  applicaple  to  all  or  most  species  are  listed 
in  this  section  (Caswell,  1989).  In  general,  reproductive  contribution  parameters  were 
calculated  as; 

Fi  =  rriiai  (2.2) 


where  rrii  is  the  maternity  function  expressed  as  female  offspring  per  female  and  ai 
is  the  probability  of  survival  to  age  one.  If  maternity  data  or  first-year  survival 
probability  were  not  available,  an  average  recruitment  of  juveniles  was  estimated 
instead  as: 

_  recruits/m^ 

adults/m2  •  ^  ’ 

In  this  case,  if  the  information  was  available,  reproductive  output  was  divided  between 
the  stages  in  proportion  to  their  contribution  to  population  reproductive  output. 
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Growth  (Gj)and  stasis  (Pj)  parameters  both  include  survival  and  both  depend  on 
stage  duration.  Growth  parameters  were  calculated  as: 

Gi  =  (2.4) 

where  7,  is  the  reciprocal  of  the  stage  duration,  and  at  is  an  estimate  of  survival 
probability  in  stage  i.  Survival  was  estimated  as 

(2.5) 

where  fx{x)  is  the  annual  mortality  rate  from  age  2;  to  age  a;  +  1  (also  called  Z^: 
in  fisheries  literature)  when  these  parameters  were  reported  in  the  literature.  The 
survival  probablity  had  to  be  estimated  in  other  ways  in  several  cases,  which  are 
detailed  with  the  parameter  estimates  below.  Stasis  parameters  were  estimated: 


Pi=(Tj(l-7i)  fori  =  1,2. 


(2.6) 


For  i  =  3, 


P,  =  exp(‘hMl)) 


where  is  the  maximum  age  reported,  and  X3  is  age  upon  entering  stage  three. 
This  function  makes  1%  of  the  population  survive  until  age  Xm-  We  identified  the 
age  which  only  1%  of  the  population  attained  from  population  age  compositions  or 
as  the  oldest  age  reported. 

Parameter  calculations  for  species  for  which  projection  matrices  were  available 
were  made  from  these  matrices.  Reported  annual  matrices  with  more  than  three 
stages  had  to  be  compressed  down  to  three  stages.  This  was  done  as  follows.  The 


proportion  of  the  population  starting  in  each  stage  was  calculated  as  the  sum  of  the 
elements  in  the  stable  stage  distribution  falling  in  the  stage.  For  each  stage,  the 
elements  outside  that  stage  in  the  stable  stage  distribution  were  replaced  with  zeros. 
This  vector  was  multiplied  by  the  annual  matrix.  Contributions  from  that  stage  to 
the  others  were  calculated  as  the  sum  of  the  elements  in  the  resulting  vector  falling 
in  each  stage,  divided  by  the  proportion  of  the  population  that  started  in  that  stage. 


44 


2.2.2  Bivalve  Species 

The  species  we  chose  to  use  were  dictated  by  the  availability  of  sufficient  information 
to  estimate  model  parameters.  Species  were  chosen  to  represent  a  range  of  average 
adult  size,  life  span,  and  habitat  latitude  and  type  (Table  2.1).  In  order  to  use  the 
same  model  structure,  which  is  useful  in  making  interspecific  comparisons  (Enright  et 
ai,  1995),  all  species  had  to  have  a  life  span  of  at  least  three  years.  Where  ages  were 
used,  we  assumed  that  shell  growth  rings  are  formed  annually  (as  was  assumed  by 
the  authors  of  literature  cited  here)  and  that  populations  were  sampled  without  bias 
with  respect  to  age.  Finally,  species  were  chosen  if  most  of  the  information  required 
to  parameterize  the  model  was  available  in  the  literature.  Where  necessary,  arbitrary 
estimates  of  missing  parameters  were  made  from  anecdotal  evidence. 

2.2.3  Matrix  Analysis 

We  calculated  population  growth  rate  as  the  dominant  eigenvalue  of  the  projection 
matrix.  Before  analyzing  them  further,  we  adjusted  all  projection  matrices  so  that 
A  was  equal  to  one.  This  was  done  by  multiplying  F2  and  F3  by  a  value,  k,  which 
changed  A  to  equal  unity.  Only  the  reproductive  parameters  were  adjusted  since  they 
are  not  only  difficult  to  measure  but  may  vary  from  year  to  year  over  wider  ranges 
than  do  growth  and  survival  (Nakaoka,  1993).  The  magnitude  of  the  required  changes 
are  also  biologically  reasonable  in  the  case  of  the  reproductive  parameters.  This  stan¬ 
dardization  is  justified  since  the  long-term  average  growth  rate  of  each  population 
must  be  close  to  one,  otherwise  species  would  tend  toward  extinction  or  overpopula¬ 
tion.  Further,  since  elasticities  change  depending  on  how  fast  a  population  is  growing 
(Brousseau  and  Baglivo,  1982),  setting  A  =  1  eliminates  differences  in  elasticity  struc¬ 
ture  due  solely  to  differences  in  population  growth  rate. 

The  remaining  analyses  were  done  on  the  adjusted  matrices.  Of  principle  interest 
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Table  2.1:  Bivalve  species  used  in  this  analysis,  listed  according  to  their  general 
life  history  traits  and  environment.  Question  marks  indicate  categories  for  which 
exmaples  could  not  be  found  in  the  literature. 


Relative  Size/Life  Span 

Large/Long _  Medium  Small/Short 


Antarctic 


none  Adamussium  colbecki  Lissarca  miliaris 

Lissarca  notorcadensis 


Temperate  Panope  abrupta 
Arctica  islandica 


Geukensia  demissa 
Yoldia  notabilis 


Gemma  gemma 
Lasaea  rubra 


Tropical  Tridacna  gigas 


? 


? 


Broadcast  Spawner 


Brooder 
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in  this  study  were  the  elasticities, 


a 


^ij  —  Sij  ^ 


ij 


(2.8) 


where 


_  dX 

^  da,, 


ViWj 


(2.9) 

<  W,  V  > 

and  V  and  w  are  the  left  and  right  eigenvectors  of  A.  Elasticities  are  proportional 
sensitivities  to  proportional  changes  in  parameter  values,  0^.  They  describe  how 
population  growth  rate  changes  with  change  in  reproduction,  growth,  and  survival, 
at  the  different  life  cycle  stages  of  a  species.  Sums  of  elasticities  to  the  Pi,  Fi,  and  Gi 
parameters  were  calculated. 

Derivation  of  the  new  method  (H.  Caswell,  pers.  comm.)  for  calculating  elastici¬ 
ties  to  (reciprocals  of  stage  durations),  nii  (reproductive  output),  and  cTj  (survival 
probabilities)  is  shown  in  the  Appendix.  These  types  of  elasticities  are  noted  as  sub¬ 
scripts  of  e;  for  example,  elasticities  to  m,  are  noted  e^j.  According  to  the  derivation, 
eigenvalue  elasticity  to  the  lower-level  maternity  parameters  is  simply  equal  to  the 
elasticity  to  the  higher-level  fertility  parameters; 

Cmi  =  e^..  (2.10) 


The  sum  of  the  elasticities  to  maternity  parameters  is  thus  equal  to  the  sum  of  the 
fertility  parameters.  The  elasticities  to  survival  sub-parameters  are  equal  to  the  sums 
of  all  parameters  in  which  the  survival  parameter  appears.  Thus, 


and  for  i  >  2, 


Cai  —  53  +  53 

j  * 


(2.11) 


e 


O'! 


(2.12) 


3 

The  sum  of  the  elasticities  of  A  to  survival  will  always  be  one,  since  each  element  of 
the  elasticity  matrix  of  A  is  included  once  in  the  calculation  of  them.  Elasticities  to 
are  calculated  as; 


—  ^Uj 


'Yij 
s  —  1 


1  dX 


(2.13) 
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where  s  is  the  number  of  stages.  This  type  of  elasticity  is  more  complex  than  the 
others  because  when  one  7  is  changed,  the  others  in  the  same  column  must  also 
change,  since  they  must  sum  to  unity.  Elasticities  to  lower-level  parameters  were  not 
calculated  for  Gemma  gemma  or  Yoldia  notabilis  because  matrices  for  these  species 
were  generated  from  matrices,  not  from  the  parameters  required  to  calculate  these 
elasticities. 

All  analyses  were  done  using  Matlab  (v.  4.2,  The  Mathworks,  Inc.). 

2.3  Parameter  estimation 

In  this  section  we  explain  how  the  data  on  each  species  were  used  to  estimate  model 
parameters.  For  each  species,  we  briefly  describe  its  natural  history,  list  the  data  that 
was  available,  and  note  how  the  stages  were  defined.  Then  we  detail  how  survival 
and  reproductive  parameters  were  calculated. 

Adamussium  colbecki-This  scallop  is  one  of  the  dominant  species  in  shallow  Antarc¬ 
tic  waters  (Stockton,  1984).  It  is  anomalous  as  an  Antarctic  species  in  its  large  size 
(12  cm)  and  unprotected,  planktotrophic  development  mode  (Berkman  et  al,  1991). 
Data  for  the  model  were  obtained  from  Berkman  (1990;  size-frequency  distributions, 
age  V.  shell  height  plot,  data  from  a  three-month  mark-recapture  survival  study) 
and  Stockton  (1984:  natural  history).  We  defined  stage  two  as  beginning  at  two 
years  of  age  because  substantial  mortality  of  small  scallops  had  occurred  between  the 
sampling  of  Stockton  (1984)  and  a  previous  population  survey  (Dayton  and  Oliver, 
1977).  We  defined  stage  three  as  beginning  at  age  five,  when  young  scallops  leave 
their  sites  of  byssal  attachment  on  the  shells  of  larger  adults  and  sexual  maturity 
is  reached  (Berkman,  1990;  S.  Van  Bloem,  1996).  Although  Berkman  (pers.  com.) 
suspects  that  scallops  live  longer  than  the  oldest  age  reported,  lacking  any  specific 
higher  value,  the  maximum  reported  age  of  20  years  (Berkman,  1990)  was  used  as  the 
estimate  of  A  survival  rate  for  stage  2  scallops  was  estimated  from  the  mortality 
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rate — three  of  twenty-five  scallops  died  in  a  three-month  period — in  a  mark-recapture 
study  (Berkman,  1990).  Assuming  a  constant  survival  rate  over  the  year,  annual  sur¬ 
vival  was  estimated  as  0.56.  Survival  in  stage  one  was  arbitrarily  estimated  as  10%  of 
that  in  stage  2,  or  0.056.  Recruitment  of  one-year  olds  was  estimated  by  first  noting 
that  one-year  old  scallops  are  approximately  15  mm  in  shell  height  (Figure  4;  Berk- 
man,  1990).  The  mean  over  two  sampling  years  of  the  ratio  of  scallops  less  than  15 
mm  to  adults  was  0.06,  according  to  size-frequency  distributions  in  the  same  paper. 
Assuming  a  1:1  sex  ratio,  recruitment  of  female  offspring  was  estimated  as  0.03. 

Arctica  islandica-The  ocean  quahog  is  an  infaunal,  subtidal  clam  inhabiting  muddy 
bottoms  in  shallow  to  deep  waters  of  the  boreal  Atlantic.  It  can  grow  to  10  cm  in 
length  and  to  an  age  of  150  years  (Thompson  et  al,  1980a).  Population  data  on 
the  ocean  quahog  were  obtained  from  Thompson  et  al.  (1980  a  and  b:  shell  length, 
shell  growth  band,  and  gonadal  condition  data),  and  Brey  et  a/. (1990:  estimates  of  Z 
calculated  as  a  regression  on  the  length-converted  catch  curve).  We  defined  stage  two 
as  starting  at  7  years  of  age  when  the  clams  reach  a  size  refuge  from  predation  and 
have  a  substantially  lower  mortality  rate  than  when  smaller  (Brey  et  al,  1990).  Stage 
three  was  defined  as  beginning  in  the  tenth  year  of  life  at  the  onset  of  sexual  maturity 
(Thompson  et  al.,  1980b).  From  Figure  4  in  Thompson  et  al.  (1980a),  Xm  was  esti¬ 
mated  as  142  years.  The  logic  of  dividing  the  population  into  these  stages  is  reinforced 
by  Brey  et  al.  (1990),  who  found  changes  in  the  slope  of  the  length-converted  catch 
curve  at  these  ages,  suggesting  that  changes  are  occurring  in  energy  allocation  and 
mortality  rates  at  these  ages.  Length-converted  catch  curves  plot  ln{Ni/At)  against 
ti,  where  Ni  is  number  of  animals  in  the  zth  length  class,  At  is  the  time  spent  in  this 
length  class,  and  U  is  age  in  length  class  i.  The  opposite  of  the  slope  of  a  regression 
fitted  through  these  points  is  an  estimate  of  total  population  mortality,  Z.  Mortality 
rates  {Zi  and  Z2)  were  obtained  from  the  legend  to  Figure  7  in  Brey  et  a/.  (1990)  and 
converted  to  survival  rates  by  Equation  2.5.  A.  islandica  recruitment  is  consistently 
described  in  the  literature  as  very  poor,  but  not  measured.  An  arbitrary  small  value. 
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0.001,  was  selected  for  F3. 

Gemma  pemma-The  gem  clam  is  a  small  (ca.  4  mm)  brooder  which  lives  in 
temperate  subtidal  sand  flats  and  has  a  life  span  of  less  than  four  years  (Weinberg, 
1985).  We  obtained  data  from  Weinberg  et  a/. (1986:  population  projection  matrices). 
All  three  stages  were  defined  as  lasting  one  year,  due  to  the  short  life  span  of  this 
clam.  All  matrix  parameters  were  calculated  from  projection  matrices  with  a  three- 
month  time  step  reported  in  Weinberg  et  al.  (1986).  The  matrix  for  1978  in  Table 
2  of  Weinberg  et  al.  (1986)  was  selected  as  representative,  and  raised  to  the  fourth 
power  to  obtain  a  projection  matrix  with  a  time  step  of  one  year.  Model  parameters 
for  this  analysis  were  calculated  from  the  one-year  matrix. 

Geukensia  demissa-The  ribbed  mussel  lives  in  the  high  intertidal  in  the  edges  of 
peat  banks,  in  salt  marshes  along  the  eastern  coast  of  the  U.S.  It  broadcast  spawns 
annually,  and  can  live  over  15  years  (Keunzler,  1961).  Data  were  obtained  from 
Keunzler  (1961;  recruits  per  adult,  size-specific  mortality  rates,  size  at  sexual  ma¬ 
turity)  and  Bertness  and  Grosholz  (1985;  size  at  age).  Sexual  maturity  is  reached 
at  about  two  years,  so  stage  two  was  defined  as  beginning  at  age  two  (Bertness  and 
Grosholz,  1985).  Stage  three  was  defined  as  beginning  at  age  four,  when  there  is  an 
increase  in  survival  rate  related  to  reaching  a  size  refuge  from  predation  (Keunzler, 
1961;  Bertness  and  Grosholz,  1985).  Fifteen  years  was  estimated  as  x^-  Survival  was 
estimated  from  Table  6  in  Keunzler  (1961)  by  first  estimating  that  mussels  in  weight 
classes  from  25-199  mg  were  in  stage  one,  and  from  200-799  mg  were  in  stage  two 
from  weight-on-size  regression  equations  (Figure  2,  Keunzler,1961).  Per  cent  mor¬ 
talities  for  the  two  model  stages  were  calculated  by  summing  the  percentages  in  each 
two-month  period  for  three  size  classes  in  each  stage,  and  dividing  by  three  hundred. 
Two-month  survival  rates  were  calculated  as  one  minus  the  mortality  rates.  Annual 
survival  for  the  two  model  stages  was  calculated  as  the  product  of  the  six  two-month 
survival  rates:  0.28  for  stage  one  and  0.71  for  stage  two.  Recruitment  was  estimated 
as  0.44,  the  sum  of  the  two-month  counts  of  recruits  per  adult,  divided  by  two  to 
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count  only  females.  It  was  partitioned  equally  between  the  two  stages  for  lack  of 
information  on  reproductive  output  by  size. 

Lasea  rubra-  This  cosmopolitan  bivalve  lives  in  the  high  rocky  intertidal  where 
it  occurs  in  Britain,  on  barnacle  shells  and  in  tufts  of  lichen.  It  is  hermaphroditic, 
and  broods  offspring  that  are  released  as  shelled  juveniles  each  year.  The  largest 
specimens  are  less  than  3  mm,  and  very  few  L.  rubra  live  past  three  years  of  age 
(Seed  and  O’Conner,  1980).  Data  were  obtained  from  Seed  and  O’Conner  (1980;  size- 
frequency  histograms,  natural  history)  and  McGrath  and  O  Foighil  (1986;  regression 
of  brood  size  on  shell  size,  reproductive  condition  data)  Stage  two  begins  at  age  two, 
and  stage  three  begins  at  age  three,  since  the  bivalve’s  life  is  too  short  to  partition  any 
differently.  Population  densities  are  known  to  fluctuate  over  short  time  scales  (Seed 
and  O’Conner,  1980),  so  estimation  of  survival  from  size-frequency  data  would  not 
be  accurate.  However,  no  other  information  was  available  on  survival,  so  arbitrary 
estimates  oi  a i  =  0.5  and  <72  =  0.8  were  made  from  observing  the  general  trend 
in  population  size-structure.  There  tend  to  be  high  numbers  of  clams  less  than  1.3 
mm,  and  smaller  but  consistent  numbers  of  larger  clams,  suggesting  that  mortality 
at  small  sizes  is  higher  than  at  large  ones.  In  order  to  estimate  fertilities,  size  ranges 
for  each  stage  were  estimated  from  size-frequency  histogram  in  Seed  and  O’Conner 
(1980).  The  May  panel  in  Figure  1  shows  cohort  sizes  at  one  year  from  the  beginning 
of  embryogenesis.  Size  ranges  were  estimated  as  0.6-1. 2  mm  for  stage  1,  1.3-1. 7  mm 
for  stage  2,  and  1.8-2. 4  for  stage  3.  According  to  Table  1  of  McGrath  and  O  Foighil 
(1986),  only  a  portion  of  the  bivalves  in  stage  2  brood.  This  portion  was  estimated 
to  be  0.26,  as  the  mean  (over  size  classes  from  1.3  to  1.8  mm)  proportion  brooding. 
Number  of  embryos  per  brood,  B,  was  estimated  for  the  median  stage  sizes,  L  =  1.5 
and  L=2.1,  from  the  regression: 

B  =  1.6L2-"«  (2.14) 

(McGrath  and  O  Foighil,  1986).  According  to  these  calculations,  stage  two  bivalves 
brood  5.2  embryos,  and  three  year  olds  brood  about  14.6  embryos  (McGrath  and 
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O’Foighil,  1986).  Multiplying  the  portion  of  stage  two  L.  rubra  brooding  by  B  yields 
1.4.  Juveniles  are  assumed  to  crawl  away  from  the  parent  on  their  first  birthday  so 
no  additional  mortality  was  imposed  on  the  brood  during  this  stage. 

Lissarca  miliaris-This  sub-littoral  Antarctic  bivalve  lives  epifaunally  on  algae  and 
other  surfaces.  It  reaches  a  maximum  size  of  6  mm,  with  a  life  span  of  up  to  7  years. 
L.  miliaris  broods  its  offspring  to  release  as  shelled  juveniles  annually.  All  of  the  data 
used  below  are  from  Richardson  (1979;  natural  history,  population  size-frequencies, 
and  size,  age,  and  reproductive  condition  data).  Embryos  are  brooded  for  an  entire 
year,  so  juveniles  enter  stage  one  just  as  they  are  released  from  the  parent.  Stage  1 
was  defined  as  the  first  year  post-brooding,  when  mortality  is  high.  Sexual  maturity 
is  not  reached  until  four  years,  so  stage  two  was  defined  as  the  pre-reproductive  period 
with  moderate  mortality-ages  2  to  4.  Stage  three  was  defined  as  beginning  at  age  four 
and  includes  all  sexually  mature  adults.  Maximum  age  was  estimated  as  six  years, 
from  the  population  size-frequency  in  Figure  12.  Richardson  (1979)  assumed  that  the 
decline  in  numbers  of  individuals  in  the  population  with  increasing  numbers  of  shell 
growth  rings  was  attributable  to  mortality.  Since  this  bivalve  is  sessile,  immigration 
and  emigration  are  minimal  and  this  assumption  is  reasonable.  The  mortality  rate 
for  stage  1  was  estimated  as  65%  in  the  first  post-brooding  year,  as  reported  by 
Richardson  (1979,  p.llO).  Richardson  (1979)  noted  that  less  than  10%  of  emerging 
juveniles  reach  age  four  (stage  3),  so  survival  from  stage  two  to  three  was  estimated 
at  0.53.  This  value  brings  the  number  of  stage  two  bivalves  to  ten  percent  of  that 
entering  stage  one  after  two  years.  Fecundity  was  estimated  by  the  mean  over  three 
months  of  the  number  of  recruits  per  adult  (1.01)  multiplied  by  0.61,  the  proportion 
of  females  in  the  population. 

Lissarca  notorcadensis-This  small,  brooding  clam  lives  byssally  attached  to  the 
spines  of  sea  urchins  and  is  abundant  on  the  Antarctic  shelf  and  slope.  Adults  reach 
a  size  of  about  8  mm.  and  live  to  about  fourteen  years  (Prezant  et  ai,  1992).  Data 
used  to  estimate  model  parameters  were  found  in  Prezant  et  al.  (1992:  natural 
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history,  reproductive  condition  and  fecundity),  Brey  and  Main  (1992:  length-at-age, 
size  frequency  distributions),  and  Brey  et  al.  (1993:  size-frequency  distributions). 
Embryos  are  brooded  for  a  year,  so  stage  one  was  defined  as  starting  at  the  time 
when  juveniles  crawl  away  from  the  parent,  making  reproductive  output  the  number 
of  shelled  embryos  brooded  to  release.  Stage  two  was  defined  as  from  age  two  to  age 
four.  Stage  three  was  defined  as  beginning  at  age  4  at  the  onset  of  sexual  maturity. 
Size  during  the  first  year  post- release  was  estimated  from  Brey  and  Hain  (Figure  3) 
as  from  1.0  mm  to  2.5  mm,  and  during  the  next  two  years,  until  sexual  maturity  is 
reached,  as  from  2.5  mm  to  4.0  mm.  The  age  beyond  which  1%  of  the  population  lives 
was  estimated  as  14  years  from  the  size-at-age  plot  in  Brey  and  Hain.  Assuming  that 
populations  have  a  stable  size  structure  and  that  there  is  no  emigration  or  immigration 
(potentially  a  reasonable  assumption  since  clams  live  byssally  attached),  decline  of 
numbers  in  progressiveley  larger  size  classes,  and  hence,  age  classes,  was  assumed 
to  be  due  to  mortality.  Size  frequency  distributions  (Brey  and  Hain  1992,  Figure  2; 
Brey  et  al,  1993,  Figure  3)  show  that  from  the  first  year  to  the  second  year,  there 
is  a  fairly  sharp  reduction  in  numbers,  but  after  that,  numbers  decline  much  more 
slowly.  In  the  population  sampled,  there  were  approximatley  14  individuals  in  the 
size  categories  falling  in  the  range  of  first  year  sizes,  and  about  4  individuals  in  the 
size  categories  for  the  next  two  year  classes,  so  survival  from  stage  1  to  stage  two  was 
estimated  as  0.29.  Survival  in  stage  2  seems  to  be  high,  since  freqencies  don’t  show 
a  decreasing  trend  during  this  time,  so  an  arbitrary  “high”  survival  rate,  0.8,  was 
chosen  for  this  value.  Fecundity  was  calculated  as  one  half  of  9.4,  the  average  brood 
size,  assuming  the  population  sex  ratio  was  1:1  (Prezant  et  al,  1992). 

Panope  abrupta-The  Pacific  geoduck  clam  is  infaunal  and  is  found  subtidally 
along  the  Pacific  coast  of  North  America  from  Baja  California  to  Alaska.  It  reaches 
sizes  up  to  20  cm  and  can  live  over  100  years.  Geoducks  broadcast  spawn  annually, 
and  population  sex  ratios  are  approximately  1:1.  Data  were  obtained  from  Sloan 
and  Robinson,  (1984;  natural  history,  histogram  of  age  frequencies,  mortality  rate) 
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and  Breen  and  Shields  (1983;  mortality  rate,  clam  densities).  Geoducks  attain  sexual 
maturity  at  about  age  five,  which  we  defined  as  the  beginning  of  stage  two.  By  age 
ten,  clams  have  essentially  stopped  growing  in  shell  length,  and  we  defined  stage  three 
as  beginning  at  that  time.  Age  beyond  which  only  1%  of  the  population  lives  was 
estimated  as  100  years  from  the  histogram  of  age-frequency  in  Sloan  and  Robinson 
(1984;  Figure  4).  Mortality  rates  were  calculated  in  these  papers  as  the  slope  of 
the  regression  of  the  natural  logarithm  of  frequency  on  age.  Mortality  rates  were 
reported  only  for  the  whole  population,  not  by  age.  Lacking  mortality  data  by  age- 
class,  we  used  the  highest  reported  mortality  rate  (0.035;  Sloan  and  Robinson,  1984) 
to  calculate  ai  (by  Equation  2.5),  and  the  lowest  reported  mortality  rate  (0.01;  Breen 
and  Shields,  1983)  to  calculate  <72.  These  values  likely  overestimate  survival  because 
only  a  small  portion  of  clams  sampled  to  generate  the  frequency  on  age  curve  were 
less  than  ten  years  old.  Breen  and  Shields  (1983)  reported  finding  35  clams  ages  1-5 
in  1,982  clams  total.  Assuming  that  those  clams  were  evenly  distributed  among  ages 
one  through  five,  there  would  be  seven  one-year  old  clams  present  in  the  sample. 
Parameters  and  Fz  were  both  estimated  as  7/1,982  =  0.0035.  There  is  likely 
to  be  a  difference  in  reproductive  output  from  clams  in  stages  two  and  three,  but 
as  there  were  no  data  addressing  this  difference,  no  attempt  was  made  to  partition 
reproductive  output  between  the  stages. 

Tridacna  pz^as-The  tropical  giant  clam,  T.  gigas  is  the  largest  and  fastest  growing 
marine  bivalve,  reaching  sizes  up  to  1  m  over  its  70  year  life  span.  They  live  only  in 
shallow  waters,  and  host  symbiotic  zooxanthellae  algae.  These  clams  are  protandric 
hermaphrodites,  and  are  broadcast  spawners  (Beckvar,  1991).  Data  were  obtained 
from  Braley  (1988;  recrutiment  rate),  Munro  (1988;  average  annual  survival  rates), 
and  Pearson  and  Munro  (1991;  average  annual  survival  rate,  length-at-age  data,  size- 
frequency  distributions).  Stage  two  was  defined  as  when  clams  become  mature  in 
the  male  phase  at  about  5  years,  and  stage  three  was  defined  as  when  clams  become 
mature  in  the  female  phase  at  about  10  years  (R.  Braley,  pers.  com.).  Maximum 
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age  reported  was  70  years  (Pearson  and  Munro,  1988;  figure  6).  Survival  rates  were 
estimated  for  sizes  of  clams  approximately  age  one  and  two  from  mean  annual  survival 
rates  (Figure  8)  in  Pearson  and  Munro  (1988)  as  0.22  and  0.31.  Survival  rates  through 
age  three  and  four  were  reported  in  Table  1  of  Munro  (1988)  as  0.505  and  0.635.  These 
rates  were  multiplied  together  to  obtain  survival  rate  through  stage  one.  Survival 
rates  from  ages  5-10  (Munro,  1988;  tablel)  were  similarly  multiplied  together  to 
obtain  survival  through  stage  two.  Natural  recruitment  in  the  field  was  measured  at 
Michaelmas  Cay,  Great  Barrier  Reef,  Australia  (Braley,  1988),  the  same  site  where 
the  survival  estimates  used  above  were  measured.  Recruitment  observed  after  one 
year  was  1.2%  of  the  population.  Only  stage  three,  the  mature  females,  reproduce  in 
this  model.  The  number  of  recruits  was  not  divided  by  two  to  account  for  proportion 
females  since  they  all  become  female  if  they  survive  long  enough. 

Yoldia  notabilis-This  clam  is  abundant  in  shallow  soft-bottom  areas  in  Japan. 
It  has  planktonic  larvae,  and  reproduces  annually  over  a  relatively  long  lifespan 
(Nakaoka,1993).  All  data  used  here  came  from  information  on  natural  history  and  a 
Leslie  matrix  for  Y.  notabilis  in  Nakaoka  (1993).  Stage  one  was  defined  as  ages  one 
through  three,  when  mortality  rates  are  high  due  to  predation.  Stage  two  was  defined 
as  beginning  at  age  four,  when  mortality  rates  are  lower  due  to  a  size  refuge  from 
predation,  and  sexual  maturity  is  reached.  We  defined  stage  three  as  beginning  at 
age  eight,  which  splits  mature  clams  into  two  groups,  the  younger  ones  having  lower 
reproductive  outputs  than  the  larger  ones.  Matrix  parameters  were  calculated  from 
the  Leslie  matrix  for  station  YA  (Table  1  in  Nakaoka,  1993),  using  the  equilibrium 
settlement  rate  for  this  site,  9.71x10“^,  as  tq. 

2.4  Results 

Estimated  parameter  values  (Table  2.2)  and  the  population  growth  rates  calculated 
from  them  (Table  2.3)  showed  that  while,  most  of  the  populations  were  not  as  near 
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equilibrium  as  we  would  expect,  poulation  growth  rates  were  at  least  reasonable. 
For  the  broadcast-spawning  bivalves,  A  was  equal  to  P3  in  some  cases  (ranging  from 
0.72  to  0.97),  while  for  the  brooding  bivalves,  A  varied  from  0.52  to  greater  than 
one  and  did  not  correspond  to  any  one  matrix  element.  For  most  of  the  broadcast 
spawning  species,  values  for  k  (Table  2.3)  were  large.  Although  population  growth 
rates  were  close  to  one,  elastisicities  (Table  2.4)  to  Fi  were  small  compared  to  those 
for  stasis  parameters,  so  large  changes  in  them  would  be  required  to  alter  A.  The  large 
elasticities  to  P3  for  broadcast-spawners  also  explained  why  this  parameter  determines 
A  for  those  species.  In  fact,  P3  always  had  the  largest  elasticity  for  broadcast  spawners 
in  this  model,  and  elasticity  to  P2  was  always  greater  than  to  Pi. 

For  any  given  species,  elasticities  to  early  and  late  growth  (Ci  and  G2)  were  nearly 
equal,  regardless  of  any  difference  in  magnitude  of  Gx  and  6*2 .  For  brooders,  elastici¬ 
ties  to  growth  and  late  reproduction  (P3)  were  the  highest,  except  for  L.  notorcadensis. 
L.  notorcadensis  lives  the  longest  of  the  brooders  analyzed — eight  years  longer  than 
the  next  most  long-lived.  Its  P3  parameter  was  larger  than  for  other  brooders,  and 
the  largest  in  the  matrix.  The  short-lived  broadcaster,  G.  demissa  showed  a  similar 
pattern,  suggesting  that  life  span,  not  reproductive  mode,  was  determining  elasticity 
structure  in  these  models. 

When  elasticities  were  summed  by  type  of  parameter,  the  largest  sum  for  all  species 
with  life  spans  of  about  14  years  or  more  (including  the  long-lived  brooder)  was  of 
stasis  parameters,  P,  while  for  the  shorter  lived  species  it  was  growth  parameters,  G. 
The  general  pattern  of  elasticity  structure  (Figure  2.4)  showed  the  smallest  brooders 
clearly  separated  from  the  largest  broadcast-spawners,  but  there  was  a  region  in  the 
center  of  the  plot  in  which  they  both  occurred. 

Calculation  of  elasticities  to  lower— level  parameters  (Table  2.4)  had  mixed  results. 
The  elasticities  to  were  the  same  as  to  the  Pj,  so  this  method  offered  no  im¬ 
provement  over  the  standard  method  for  those  parameters.  Since  the  elasticities  to 
cr,  always  sum  to  one,  comparing  their  sum  between  species  is  useless.  However, 
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Table  2.2:  Parameter  Estimates.  Estimated  parameter  values  for  calculation  of  matrix 
elements:  and  2:3  are  ages  at  the  beginning  of  stages  two  and  three,  and  is 

the  maximum  age,  7i  are  reciprocals  of  stage  durations,  and  Oi  are  survival  rates; 
and  calculated  matrix  elements:  Gj,  and  Pj.  Dashes  indicate  parameters  which 

were  not  estimated  for  a  species.  The  model  for  G.  gemma  includes  two  additional 
parameters:  Gi,3  =  0.03,  and  P2,3  =  0.452. 


Species 

X2 

2:3 

Ti 

72 

0-1 

0'2 

A.  colbecki 

2 

6 

20 

1.0 

0.25 

0.056 

0.56 

A.  islandica 

7 

10 

142 

0.143 

0.333 

0.633 

0.942 

G.  gemma 

2 

3 

4 

- 

- 

- 

- 

G.  demissa 

2 

4 

15 

1.0 

0.5 

0.28 

0.71 

L.  rubra 

2 

3 

4 

1.0 

1.0 

0.5 

0.8 

L.  miliaris 

2 

4 

6 

1.0 

0.5 

0.35 

0.53 

L.  notorcadensis 

2 

4 

14 

1.0 

0.5 

0.29 

0.8 

P.  abrupta 

5 

10 

100 

0.2 

0.2 

0.97 

0.99 

T.  gigas 

5 

10 

70 

0.2 

0.0167 

0.0219 

0.0465 

Y.  notabilis 

4 

8 

15 

- 

- 

- 

- 

Species 

Pi 

P2 

P3 

Gi 

G2 

F2 

P3 

A.  colbecki 

0 

0.375 

0.72 

0.056 

0.185 

0 

0.03 

A.  islandica 

0.542 

0.628 

0.966 

0.091 

0.314 

0 

.001 

G.  gemma 

0.296 

0.134 

0 

0.09 

00035 

1.62 

0.452 

G.  demissa 

0 

0.355 

0.66 

0.28 

0.355 

0.22 

0.22 

L.  rubra 

0 

0 

0.01 

0.5 

0.8 

1.4 

14.6 

L.  miliaris 

0 

0.265 

0.1 

0.35 

0.265 

0 

0.62 

L.  notorcadensis 

0 

0.4 

0.63 

0.29 

0.4 

0 

4.7 

P.  abrupta 

0.776 

0.792 

0.95 

0.194 

0.198 

0.0035 

0.0035 

T.  gigas 

0.0175 

0.456 

0.926 

0.00438 

0.00777 

0 

0.012 

Y.  notabilis 

0.07 

0.68 

0.0797 

0.003 

0.17 

33.12 

7.25 
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Table  2.3:  Annual  population  growth  rates,  A,  of  bivalve  populations  (before  adjusting 
parameters),  and  constant,  k,  used  to  adjust  fertility  parameters  prior  to  calculating 
elasticities.  _ 


Species 

A 

k 

A.  colbecki 

0.72 

563 

A.  islandica 

0.97 

204 

G.  demissa 

0.75 

2.88 

G.  gemma 

0.66 

3.38 

L.  rubra 

1.93 

0.15 

L.miliaris 

0.52 

11.5 

L.  notorcadensis 

1.19 

0.4 

P.  abrupta 

0.95 

13.8 

T.  gigas 

0.92 

9.5x10' 

Y.  notabilis 

0.82 

2.88 
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Table  2.4:  Eigenvalue  elasticities  of  adjusted  matrices.  Matrix  elements  are  found 
in  Table  2,  but  Fi  have  been  multiplied  by  k  (Table  3)  prior  to  calculating  elastici¬ 
ties.  For  placement  of  each  element  in  the  displayed  arrays,  see  definition  of  A  with 
Equation  2.1. 


A.  colbecki 

0 

0 

0.1620 

0.1620 

0.0972 

0 

0 

0.1620 

0.4167 

A.  islandica 

0.0345 

0 

0.0292 

0.0292 

0.0486 

0 

0 

0.0292 

0.8293 

G.  demissa 

0 

0.1207 

0.1260 

0.2467 

0.1358 

0 

0 

0.126 

0.2447 

G.  gemma 

0.1562 

0.3469 

0.0246 

0.3018 

0.0538 

0.0461 

0.0679 

0.001 

0 

L.  rubra 

0 

0.0365 

0.3079 

0.3445 

0 

0 

0 

0.3079 

0.0031 

L.  miliaris 

0 

0 

0.2880 

0.2880 

0.1039 

0 

0 

0.2880 

0.0320 

L.  notorcadensis 

0 

0 

0.1862 

0.1862 

0.1242 

0 

0 

0.1862 

0.3171 

P.  abrupta 

0.1373 

0.0080 

0.0316 

0.0396 

0.1509 

0 

0 

0.0316 

0.6010 

T.  gigas 

0.0011 

0 

0.0610 

0.0610 

0.0530 

0 

0 

0.0610 

0.7630 

Y.  notabilis 

0.0177 

0.2265 

0.0092 

0.357 

0.5009 

0 

0 

0.0092 

0.0008 
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Figure  2.4;  Elasticity  Structures  of  Bivalve  Life  Histories.  The  sums  of  elasticities  to 
Pii  P'iy  &rid  Gi  plotted  against  one  another.  The  squares  are  brooding  species;  the 
circles  are  broadcasting  species.  Note  that  only  two  of  the  variables  are  independent 
as  all  elasticities  sum  to  one. 
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Table  2.5:  Elasticities  to  lower-level  parameters.  For  each  species,  elasticities  to  stage 
transition  probability,  7ij,  and  survival  probability,  ai  are  shown.  Elasticities  to  rrii 
are  not  shown  because  they  are  the  same  as  elasticites  to  Fi  (Table  2.4). 


Species 

Ol 

0-2 

(^3 

A.  colbecki 

0.4792 

0.3241 

0.2593 

0.4167 

A.  islandica 

0.3716 

0.0929 

0.0778 

0.8293 

G.  demissa 

0.3415 

0.4935 

0.2618 

0.2447 

L.  rubra 

0.2378 

0.689 

0.3079 

0.00511 

L.  miliaris 

0.0604 

0.5761 

0.3919 

0.032 

L.  notorcadensis 

0.3862 

0.3725 

0.3104 

0.3171 

P.  abrupta 

0.3405 

0.2165 

0.1825 

0.0601 

T.  gigas 

-9.9104 

0.123 

0.114 

0.763 

comparing  the  pattern  of  ai  values  from  one  species  to  another  had  similar  results 
to  the  standard  elasticity  analysis:  elasticity  to  ai  was  highest  for  brooders,  while 
for  most  broadcasters,  elasticity  to  <73  was  highest.  As  we  calculated  them,  it  was 
possible  for  elasticities  to  7j  to  be  negative,  indicating  that  elasticity  to  other  param¬ 
eters  in  the  same  column  is  higher.  They  also  measured  elasticity  to  stage  durations, 
which  in  these  models  were  set  rather  artificially.  So  in  this  example,  elasticites  to  7 
parameters  were  difficult  to  interpret. 


2.5  Discussion 

There  is  a  continuum  of  contributions  to  A  by  stasis  parameters  for  the  species  ana¬ 
lyzed  here,  but  as  strategy  changes  from  brooding  and  broadcast  spawning  the  pattern 
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shifts  from  fitness  depending  on  growth  and  fecundity  to  it  depending  on  stasis.  Our 
results  suggest  that  the  longer  the  life  span,  the  less  important  reproduction  and  early 
growth  are  in  the  life  history,  regardless  of  reproductive  mode.  Menge  (1977)  made  a 
back-of-the-envelope  calculation  to  show  that  the  small,  brooding,  asteroid,  Leptas- 
terias  hexactis  would  have  to  have  a  life  span  of  over  1500  years  to  be  successful  as 
a  broadcast  spawning  species.  This  is  due  to  the  small  volume  available  for  egg  pro¬ 
duction  and  low  survival  of  planktonic  larvae.  While  it  makes  sense  for  small  size  to 
constrain  bivalves  to  brooding,  it  is  less  clear  why  larger  sized  animals  do  not  brood. 
Life  history  analyses  on  ophiuroids  have  also  shown  that  even  for  a  brooding  species, 
selection  pressure  appears  to  be  towards  long  adult  life  span  (Medeiros-Bergen  and 
Ebert,  in  prep).  After  comparing  various  asteroid  life  histories,  Ebert  (1996)  sug¬ 
gested  that  an  important  factor  in  reproductive  mode  is  the  negative  relationship 
between  larval  survival  and  life  span. 

Although  parameter  estimates  made  for  this  analysis  were  from  data  that  were,  in 
most  cases,  not  collected  for  demographic  modeling,  it  was  still  possible  to  calculate 
approximate  matrix  parameters.  A  lack  of  appropriate  data  is  blamed  for  the  limited 
analyses  that  have  been  done  to  compare  life  history  strategies,  but  we  have  found  that 
although  the  data  were  not  ideal,  there  are  more  available  than  has  been  previously 
utilized.  We  found  that  while  fecundity  data  were  commonly  reported,  growth  and 
survival  were  not  generally  measured  for  different  size  classes,  and  that  the  biggest 
gap  in  the  literature  was  in  information  on  juvenile  survival  and  recruitment.  The 
results  presented  here  offer  an  improvement  over  the  purely  statistical  comparisons 
of  life  history  traits  as  done,  for  example,  by  Way  (1988),  in  that  they  combine  the 
different  elements  of  a  life  history  strategy  and  examine  them  as  an  integrated  unit, 
rather  than  as  separate  parts. 

Growth  rates  of  bivalves  are  known  to  vary  depending  on  their  microhabitat  {e.g. 
Bertness  and  Grosholz,  1985).  Survival  {e.g.  Weinberg,  1985)  and  recruitment  {e.g. 
Braley,  1988)  are  also  known  to  vary  spatially  and  temporally.  The  projection  ma- 
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trices  calculated  here  are  mere  snapshots  of  a  range  of  possible  population  processes 
for  these  species.  Our  object  in  calculating  them  was  only  to  capture  the  salient  life 
history  features  in  general.  Studies  that  have  examined  the  effects  of  variation  in 
matrix  elements  have  shown  that  including  variation  can  lead  to  changes  in  popula¬ 
tion  growth  rates  (reviewed  in  Nakaoka,  1996).  However,  if  we  assume  that  variation 
acts  similarly  on  the  different  species  compared  here,  it  can  be  ignored  and  the  life 
histories  can  be  compared  as  outlines  representative  of  types  of  species. 

Enright  et  al.  (1995)  demonstrated  that  either  holding  the  number  of  stages 
constant  between  species  or  holding  the  time  spent  in  stages  constant  between  species 
allows  meaningful  comparisons  to  be  made  between  species.  Limited  data  forced  us  to 
hold  the  number  of  stages  constant,  and  this  approach  tends  to  increase  the  apparent 
relative  importance  of  stasis  in  the  life  history  because  more  individuals  remain  in 
each  stage  every  time  step,  and  less  proceed  to  the  next  stage  (Enright  et  al,  1995). 
Our  analysis  does  pinpoint  life  span  as  a  potentially  important  life  history  trait. 
However,  the  alteration  in  elasticities  due  to  collapsing  life  cycles  to  fewer  stages  is 
not  large  enough(see  Figure  3  in  Enright  et  al,  1995)  to  alter  the  general  patterns, 
and  it  tends  to  shift  all  species  in  the  same  direction. 

Elasticity  patterns  fall  along  a  line  over  only  a  portion  of  the  range  in  the  sum  of 
elasticity  to  E).  One  might  have  expected  that  brooding  species  would  have  higher 
elasticities  to  reproduction,  since  juvenile  survival  is  higher,  but  this  was  not  the 
case.  Reproductive  output  did  not  contribute  the  most  to  fitness  in  either  type  of 
life  history.  It  is  possible  that  bivalves  are  architecturally  constrained  in  their  life 
history  strategies  to  this  region  of  the  plot.  However,  we  only  used  ten  species  in 
this  analysis-a  minute  sample  size.  The  exclusion  of  species  with  life  spans  shorter 
than  three  years  definitely  biased  the  outcome  away  from  the  bottom  left  corner  of 
Figure  2.4,  where  reproduction  is  the  most  important  factor.  More  data,  especially 
on  tropical  species,  which  can  grow  very  fast  but  have  short  life  spans,  would  be  likely 
to  fill  in  the  gaps  of  diversity  of  bivalve  life  histories. 
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Analyses  on  plants  (Silvertown  et  ai,  1993)  have  shown  that  they  cover  more 
regions  of  the  elasticity  space  than  do  bivalves,  however,  the  plants  analyzed  represent 
many  classes  and  the  species  analyzed  here  are  all  in  Bivalvia.  The  plants  clearly  show 
distinct  ecological  groupings  in  demographic  G-P-F  space,  in  a  similar  manner  to  the 
patterns  demonstrated  by  Grime  (1977)  in  trait  G-S-R  space.  Brooding  bivalves  fall 
into  the  same  region  in  elasticity  space  as  do  semelparous  and  iteroparous  herbs,  while 
the  larger  broadcasters  fall  into  the  same  region  as  do  woody  plants  (Silvertown  et 
al,  1993).  Silvertown  suggests  that  the  G-P-F  method  is  valuable  because  fitness 
consequences  of  the  life  history  can  be  directly  interpreted  from  it. 

The  method  described  here  to  calculate  elasticites  to  lower-level  parameters  had 
mixed  results.  We  think  that  elasticities  to  oi  are  potentially  useful  in  comparative 
demography,  while  elasticities  to  rrii  and  7^-  were  not  informative.  Because  of  dis¬ 
satisfaction  with  the  fact  that  elasticites  to  P  and  G  parameters  confound  survival 
and  growth,  it  is  worth  pursuing  alternative  methods  of  calculating  elasticities  to 
lower-level  parameters  further. 

One  of  the  main  points  of  Grime’s  triangle  scheme  is  that  the  three  types  of  life 
history  strategies  are  found  in  species  occupying  particular  habitats.  The  link  between 
the  environment  and  the  strategy  of  the  plant  is  assumed  to  be  tight.  For  bivalves, 
this  connection  is  not  obvious.  Thorson  (1950)  suggested  that  environmental  factors 
at  different  latitudes  and  water  depths  determine  the  reproductive  strategy  of  marine 
invertebrates,  but  many  counterexamples  have  been  discovered  more  recently.  The 
proportion  of  bivalve  species  that  brood  may  be  higher  in  communities  in  the  Antartic 
and  the  deep  sea  (Grahme  and  Branch,  1985),  but  other  life  histories  are  present  as 
well.  Brooding  and  broadcasting  species  are  found  in  all  environments,  whether 
categorized  by  latitude  or  by  tidal  level.  The  tropical  giant  clam  has  essentially 
the  same  contributions  to  fitness  from  different  parts  of  its  life  history  strategy  as 
the  temperate  geoduck  clam  or  the  Antarctic  scallop.  Brooders  are  found  in  the 
high  intertidal,  a  very  stressful  environment,  but  so  are  broadcasters.  We  do  not 
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think  that  the  aspect  of  comparative  life  history  strategy  dealing  with  correlations  to 
environmental  types  is  as  applicable  to  marine  species  as  to  plants. 

Key  differences  in  brooding  and  broadcast  spawning  life  histories  may  relate  to 
time  scales  of  variability.  The  longer  the  life  span,  the  less  important  reproductive 
success  is  in  each  year.  Long  life  spans  may  have  led  to  the  possibility  of  lower 
investment  into  riskier  reproductive  patterns.  On  the  other  hand,  reduced  energy 
investment  in  reproduction  may  allow  long  life.  Other  authors  have  contemplated 
this  marked  difference  in  reproductive  mode  by  asking  the  questions:  “what  are  the 
costs  and  benefits  of  having  planktonic  larvae?”  (Grahme  and  Branch,  1985),  “brood 
or  broadcast?”  (Menge,  1975),  or  “to  brood  or  not  to  brood?”  (Strathmann  et 
al,  1984).  Perhaps  the  question  to  ask  is:  “How  long  to  live?”  The  answer  to  this 
quesion  may  depend  on  the  regime  of  variability  in  the  environment  (Strathmann  and 
Strathmann,  1982).  Amount  of  selection  pressure  on  life  span  may  also  contribute 
to  evolutionary  outcomes.  Further  work  incorporating  variability  into  demographic 
models  on  life  history  strategies  may  be  able  to  address  this  question. 


2.6  Summary 

1.  We  suggest  that  there  are  more  data  available  for  comparative  demography  than 
have  been  utilized,  and  that  although  the  data  may  not  have  been  collected  for 
parameterizing  models,  they  are  adequate  for  theoretical  analyses. 

2.  Survival  is  important  to  both  brooders  and  broadcasters,  but  long  life  span  con¬ 
tributes  the  most  to  fitness  in  broadcast  spawning  bivalves,  and  growth  is  more  im¬ 
portant  to  the  brooders.  There  is  a  region  of  elasicity  space  where  both  brooders  and 
broadcasters  are  found. 

3.  Patterns  of  elasticity  distributions  in  marine  bivalves  apparently  do  not  correlate 
to  environmental  conditions  as  do  those  of  plants. 


65 


4.  A  new  method  for  calculating  separate  elasticities  to  growth  and  survival  was 
derived  and  used  on  these  models. 
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2.8  Appendix 

Derivations  of  elasticities  to  the  lower-level  parameters  mi,  ^ij,  and  (from  which 
projection  matrix  parameters  Pi,  Gi,  and  Fi  are  calculated)  are  shown  below  (H. 
Caswell,  pers.  comm.).  These  elasticities  were  derived  in  order  to  try  to  calculate 
elasticities  to  survival  separately  from  growth  and  reproduction,  which  is  not  possible 
by  calculating  elasticities  to  projection  matrix  parameters. 

Let 

A  =  T-hF 
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For  (jj  where  2  =  1,  sensitivities  are 


dX  ^  dX  dFj  ^  dX  dtji 

dai  dFj  day  '  dtji  day 

(2.17) 

^  dX  ^  dX 

(2.18) 

For  2  >  2, 

dX  dX  dtji  dX 

dai  j  dtji  dai  j  dtjP^'' 

(2.19) 

Sensitivities  for  are 

simply 

Sensitivities  to  jij  are: 
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(2.20) 
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Elasticities  to  different  parameters  are  noted  as  subscripts  to  the  letter  e,  for  example 


Crr,,  are  elasticities  to  m,-. 


For  i  =  li  elasticities  to  are 
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For  i  >2,  elasticities  to  ai  are 


(2.23) 
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For  'jij,  elasticities  calculated  in  the  same  manner  as  those  above  are; 
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However,  when  one  7  is  changed,  the  others  in  the  same  column  must  also  change, 
since  they  sum  to  one,  so  Equation  2.25  is  not  entirely  correct.  Suppose 


^Jkj  = 


^ij 
S  —  1 


(2.26) 


where  s  is  the  number  of  stages.  Then 
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The  role  of  recruitment  variability 
in  determining  population  growth 
rate  of  the  soft-shell  clam 

3.1  Introduction 

The  life  history  of  the  clam  Mya  arenaria,  like  that  of  many  benthic  marine  inverte¬ 
brates,  includes  a  sessile  adult  stage  that  broadcasts  large  numbers  of  gametes  into 
the  water  column.  A  female  M.  arenaria  may  shed  up  to  10®  eggs  in  the  spring  and 
the  fall  (Brousseau,  1978),  but  high  mortality  rates  and  transport  away  from  the  site 
can  result  in  undetected  settlement  of  spat,  despite  such  high  reproductive  output 
(Brousseau  et  al,  1982).  Stochasticity  in  environmental  factors,  such  as  wind  forcing 
of  currents,  that  affect  planktonic  larvae  can  cause  dramatic  spatial  and  temporal 
variability  in  settlement  {e.g.  Goshima,  1982;  Vahl,  1982;  Bachelet,  1986).  Mortality 
rates  of  planktonic  larvae  also  vary  interannually  due  to  their  dependence  on  multiple 
biological  and  physical  factors  (Rumrill,  1990). 

Variabilty  in  recruitment  is  an  important  phenomenon  in  and  of  itself,  rather  than 
just  an  irritant  to  field  biologists  (Fogarty,  1993a).  We  loosly  define  “recruits”  as  re- 


74 


cently  settled  or  obviously  young-of-the-year  individuals.  Uncertain  juvenile  survival 
is  generally  found  in  organisms  with  iteroparity  and  high  fecundity  (Murphy,  1968). 
M.  arenaria  lives  for  up  to  twelve  years  (Newell,  1991),  so  successful  recruitment  is 
not  necessary  each  year  to  maintain  populations  (Goodman,  1984).  This  combination 
of  life  history  traits  results  in  what  is  termed  the  “storage  effect,”  where  population 
growth  rate  is  more  strongly  affected  by  conditions  favorable  to  recruitment  success 
than  by  unfavorable  ones.  This  effect  has  been  shown  to  be  important  in  the  context 
of  interspecific  competition  {e.g.,  Warner  and  Chesson,  1985).  In  this  paper,  we  ex¬ 
amine  the  effects  of  recruitment  variability  on  population  dynamics  of  M.  arenaria, 
using  a  stochastic  matrix  population  model. 

Matrix  population  models  have  been  applied  to  M.  arenaria  to  calculate  sensi¬ 
tivity  to  changes  in  vital  rates  (Brousseau  and  Baglivo,  1984),  and  to  assess  disease 
impact  (Weinberg  et  al,  1996)  and  fisheries  management  strategies  (Malinowski  and 
Whitlatch,  1988).  These  studies  have  all  used  deterministic  models  that  neglect  the 
effects  of  variability.  Brousseau  and  Baglivo  (1984)  and  Malinowski  and  Whitlatch 
(1988)  used  an  assumption  of  constant  population  size  to  estimate  the  “equilibrium 
settlement  rate.”  This  value  is  defined  as  the  proportion  of  eggs  spawned  that  must 
successfully  develop  and  survive  to  their  first  year  to  produce  a  constant  population 
size  when  combined  with  all  the  other  vital  rates.  Settlement  in  this  sense  is  what 
we  examined  to  study  variation  in  recrutiment.  Since  the  population  is  assumed 
to  be  at  equilibrium  in  these  models,  only  equilibrium  settlement  rate  can  be  com¬ 
pared  among  sites,  conditions,  or  studies,  not  population  growth  rates.  Brousseau  et 
aZ.(1982)  calculated  that,  for  a  Cape  Anne,  Massachusetts  population  of  M.  arenaria, 
only  0.0001%  of  eggs,  or  about  one  recruit  or  settler  for  every  ten  adults,  needed  to 
survive  each  year  to  maintain  constant  population  size.  No  set  at  all  was  detected  in 
two  out  of  the  three  years,  however,  so  in  order  to  maintain  the  population,  a  much 
higher  settlement  must  occur  occasionally. 

To  compare  previous  work  to  the  stochastic  model  presented  here,  we  first  formu- 
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lated  a  deterministic  model.  A  deterministic  matrix  population  model  can  be  written 
as 

n(t  +  l)  =  An(t)  (3.1) 

where  n  is  a  vector  of  stage  abundances  and  A  is  a  constant  projection  matrix.  The 
projection  matrix  elements  are  estimates  of  the  possible  transitions  between  stages, 
over  one  time  step.  The  time  step  chosen  is  often  one  year.  An  annual  projection 
matrix  can  also  be  constructed  from  matrices  describing  seasonal  transitions.  For 
example,  a  population  can  be  projected  from  one  spring  to  the  next  by  a  matrix: 

spring -^winter Jail -^summer  (3-2) 

where  the  seasonal  matrices  describe  vital  rates  over  portions  of  the  year.  For  example, 
■^spring  dcscribes  growth  and  survival  from  March  to  June.  The  other  matrices  are 
similarly  defined  over  three  month  intervals.  This  approach  makes  it  possible  to 
include  data  on  a  finer  than  annual  time  scale  into  a  model  with  an  annual  time  step. 

The  dynamics  of  the  population  are  determined  by  the  matrix  A  (Caswell,  1989). 
The  dominant  eigenvalue  of  A,  A,  gives  the  population  growth  rate.  The  correspond¬ 
ing  right  eigenvector  w  and  left  eigenvector  v,  defined  by 

Aw  =  Aw  (3.3) 

v'A  =  Av'  (3.4) 

give  the  stable  stage  distribution  and  reproductive  value  distributions.  It  is  sometimes 
convenient  to  describe  population  growth  by  the  continuous  time  growth  rate.  In  A, 
also  known  as  r,  the  instantaneous  rate  of  increase. 

The  sensitivity  of  A  to  changes  in  the  entries  of  A  is  given  by 

dX/daij  =  ViWj/  <  V,  w  >  (3.5) 

where  <  v,w  >  is  the  scalar  product  of  v  and  w.  The  proportional  sensitivity,  or 
elasticity,  of  A  is  given  by 

Gij  ctijSij  f  X.  (3.6) 
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In  a  stochastic  matrix  model,  some  elements  of  A  vary  randomly,  so  that  the 
model  is: 

n(t  + 1)  =  Atn(i).  (3.7) 

Specifying  such  a  model  requires  a  stochastic  model  for  the  elements  of  A^;  what  their 
probability  distributions  are,  whether  they  are  correlated  or  vary  independently,  and 
how  the  values  at  one  time  are  related  to  those  at  proceeding  times.  Although  a  de¬ 
tailed  stochastic  model  could  require  much  more  data  than  a  deterministic  one,  simple 
approximations  of  stochastic  processes  can  be  used  to  simulate  natural  variability  for 
a  variety  of  demographic  factors  (Tuljapurkar  and  Caswell,  1996). 

There  are  now  powerful  methods  available  to  analyze  stochastic  models  in  an 
analagous  fashion  to  more  familiar  analyses  of  deterministic  models(Caswell,  1989; 
Tuljapurkar  1990;  Tuljapurkar  and  Caswell,  1997).  Population  growth  is  described 
by  the  stochastic  growth  rate  A^.  It  is  not  calculated  as  an  eigenvalue,  but  as  a 
long-term  average  of  the  growth  rates  along  a  trajectory  produced  by  a  simulation 
of  Equation  3.7.  The  stochastic  growth  rate  describes  the  long-term  average  of  every 
such  realization  with  probability  one,  and  measures  fitness  in  stochastic  life  history 
models  (Tuljapurkar  1990).  Sensitivities  of  In  to  changes  in  matrix  elements  can 
also  be  calculated  for  the  stochastic  model,  using  the  stochastic  analogs  of  the  stable 
stage  distribution  and  reproductive  value  vectors. 

Examples  of  studies  using  stochastic  models  include  those  randomly  selecting  one 
of  several  possible  matrices  representing  differen  environmental  conditions  (Silva  et 
ah,  1991;  Aberg,  1992;  Noda  and  Nakao,  1996),  those  generating  a  series  of  random 
matrices  by  jack-knifing  the  set  of  measured  vital  rates  (Gotelli,  1991),  and  those  se¬ 
lecting  a  single  parameter  randomly  from  a  specified  distribution  (Cohen  et  al,  1983; 
Nakaoka,  1997).  The  last  approach  has  the  advantage  of  generating  stochasticity  in 
the  model  from  a  continuous  distibution,  which  includes  a  continuous  range  of  values 
rather  than  just  a  few  points,  and  is  the  approach  we  chose  to  take. 

For  M.  arenaria,  the  variability  in  larval  settlement  and  survival  to  age  one  is 
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vast  compared  to  variability  in  other  matrix  elements  (Goshima,  1982).  Nakaoka 
(1997)  showed  that  incorporating  stochasticity  in  recruitment  into  a  matrix  model 
had  a  greater  impact  on  than  did  incorporating  stochasticity  in  shell  growth  rates. 
Eigenvalue  sensitivity  to  settlement  rate  has  also  been  shown  to  be  higher  than  sen¬ 
sitivity  to  adult  survival  in  M.  arenaria  (Brousseau  and  Baglivo,  1984).  For  these 
reasons,  and  because  we  were  particularly  interested  in  the  process  of  settlement,  we 
chose  to  vary  only  the  settlement  parameter  in  this  study.  Settlement  variability  was 
modelled  according  to  a  lognormal  distribution,  which  mimics  the  high  degree  of  vari¬ 
ability  in  settlement  of  benthic  marine  invertebrates  that  is  commonly  observed  {e.g. 
Seed  and  Brown,  1975;  Gaines  and  Bertness,  1993)  better  than  a  normal  distribution. 

The  goal  of  this  paper  is  to  determine  the  effect  of  stochastic  settlement  in  a 
matrix  population  model,  where  the  settlement  is  lognormally  distributed.  We  wanted 
to  compare  the  stochastic  model  to  the  results  of  the  deterministic  case,  which  has 
been  previously  analyzed  for  marine  invertebrates  {e.g.,  Malinowski  and  Whitlatch, 
1988).  This  analysis  extends  the  work  of  Nakaoka  (1997)  on  the  marine  bivalve  Yoldia 
notabilis,  who  showed  that  the  type  of  recruitment  distribution  creating  stochasticity 
in  a  model,  and  the  choice  of  parameters  to  vary  greatly  influence  the  population 
growth  rate  calculated. 

3.2  Methods 

3.2.1  Overview 

We  measured  growth  and  survival  of  adult  M.  arenaria  and  gathered  other  informa¬ 
tion  required  for  the  model,  such  as  population  size  structure,  in  Barnstable  Harbor, 
Massachusetts.  These  data  provided  all  of  the  information  to  parameterize  the  model 
except  for  reproductive  contributions.  Settlement  for  both  the  deterministic  and 
stochastic  models  was  estimated  from  data  on  M.  arenaria  populations  in  Japan 
(Goshima,  1982),  because  a  larger  set  of  data  on  spatial  and  temporal  variation  in 
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settlement  was  available  in  this  paper  than  from  our  field  data.  While  this  combi¬ 
nation  of  data  from  different  sites  does  not  allow  us  to  describe  the  states  of  the 
clam  populations  in  Japan,  Massachusetts,  or  anywhere  else,  it  allows  us  to  explore 
theoretically  the  importance  of  varying  settlement  against  the  background  of  a  set 
of  reasonable  vital  rates.  We  assume  that  these  parameters  are  accurate  enough  to 
describe  in  general  the  population  dynamics  of  M.  arenaria. 

To  investigate  the  role  of  variability  in  settlement,  we  formulated  a  stochastic 
model.  We  chose  to  include  stochastic  variation  in  only  larval  settlement,  holding 
the  rest  of  the  matrix  fixed.  The  analyses  cannot  actually  distinguish  whether  the 
variation  is  caused  by  changes  in  fertility  ,  settlement,  or  larval  survival  since  all  are 
included  in  calculation  of  these  parameters,  but  we  chose  to  interpret  the  results  as 
if  all  the  variation  is  due  to  changes  in  settlement.  In  the  absence  of  information  to 
specify  a  more  complex  temporal  pattern,  we  modelled  settlement  as  varying  indepen¬ 
dently  from  year  to  year.  This  means  that  we  needed  to  specify  only  the  distribution 
from  which  reproductive  parameter  values  were  drawn.  The  deterministic  case  was 
modeled  as  a  special  case  of  the  stochastic  one,  in  which  the  variance  of  the  settlement 
distribution  was  zero. 

3.2.2  Estimation  of  Adult  Growth  and  Survival 

Growth  and  survival  probabilities  were  measured  in  a  mark-recapture  study  in  Barn¬ 
stable  Harbor,  Massachusetts  (similar  to  Weinberg  et  ai,  1996).  In  March,  1995, 
sufficient  numbers  of  clams  were  dug  from  the  intertidal  mudfiat  to  select  twenty 
in  each  of  the  following  shell  length  size  classes;  20-40  mm  (class  1)  ;  40-49.9  mm 
(class  2)  ;  50-59.9  mm  (class  3)  ;  60-69.9  mm  (class  4)  ;  and  >70  mm  (class  5).  Size 
classes  were  selected  to  evenly  span  the  sizes  present  in  the  population.  Although 
clam  growth  slows  with  size,  enough  large  clams  grew  fast  enough  that  we  could  mea¬ 
sure  transitions  between  the  larger  size  classes.  Each  clam  was  measured  and  marked 
with  an  identification  number  on  both  valves  in  indelible  pen.  Ten  clams  (two  from 
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each  size  class)  were  placed  together  in  each  of  ten  sediment-filled  mesh  bags  (mesh 
size  ca.  20  mm),  in  their  normal  orientation  with  siphon  up  and  at  the  appropriate 
depth  for  their  size  (about  one  body  length).  The  bags  did  not  exclude  predators, 
but  simply  marked  the  location  of  the  clams  and  allow  for  easier  sampling.  The 
density  of  clams  within  the  bags  (80/m^)was  not  greater  than  natural  densities  (ca. 
400/m^).  In  June,  September,  and  December,  the  bags  were  removed,  and  the  clams 
were  sorted  from  the  sediment,  counted,  and  measured.  At  each  sampling  time,  new 
clams  were  collected  to  replace  missing  or  dead  clams,  labelled,  and  mixed  with  the 
survivors  before  placing  all  clams  in  new  bags  as  in  March  1995.  Surviving  clams  that 
were  re-deployed  were  assumed  to  have  the  same  survival  and  growth  rates  as  newly 
deployed  clams.  In  March  1996,  all  clams  were  removed  from  the  bags,  counted,  and 
measured. 

3.2.3  Model  Formulation 

The  life  cycle  graph  (Figure  3.1),  and  the  corresponding  projection  matrix  A  (Fig¬ 
ure  3.2),  depict  the  structure  of  the  model  and  its  parameters.  The  model  is  stage 
classified,  with  birth-pulse  reproduction  and  a  pre-breeding  census.  Stages  1-5  are 
the  size-classes  for  which  growth  and  survival  were  measured  in  the  field.  Since  adult 
growth  and  survival  data  were  collected  over  four  contiguous  three-month  periods, 
parameters  for  seasonal  projection  matrices  were  calculated.  First,  transition  tables 
(T)  were  constructed  from  these  data  (Table  3.1),  then  the  elements  of  the  seasonal 
projection  matrices,  aij,  were  calculated  from  the  T: 

dij  ^ij  /  “dj  f  (3.8) 

(Caswell,  1989),  where  i  is  the  row  and  j  the  column  of  the  element  in  A  or  T,  and 
rij  is  the  number  of  clams  deployed  minus  the  number  of  missing  clams.  The  seasonal 
matrices  (Table  3.2)  were  multiplied  together  to  make  a  single  projection  matrix  with 
an  interval  of  one  year  according  to  Equation  3.2. 
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Figure  3.1:  The  life  cycle  graph  for  M.  arenaria.  The  circles  indicate  stages,  and 
the  arrows  represent  these  transitions:  surviving  and  growing  to  a  subsequent  stage 
{Gij),  surviving  and  staying  in  the  same  stage  {Pi),  or  reproductive  contribution  to 
the  first  stage  (Ri).  These  parameters  are  denoted  Ri,  rather  than  Fi,  because  they 
measure  /Recruitment  at  the  end  of  the  model  time  step,  rather  than  Fertility  at  the 
beginning  of  the  time  step.  Stages  A,  B,  and  C  are  dummy  stages,  as  defined  in 
Model  Formulation. 
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Figure  3.2:  Elements  of  the  seasonal  projection  matrices,  A^,  corresponding  to  the 
stage  transitions  in  Figure  3.1.  All  locations  in  the  matrix  where  no  parameter  is 
shown  are  zero.  In  this  model,  Ga,Bi  Gb,c,  and  Gc,i  are  always  one  and  are  not 
displayed  in  subsequent  tables  of  seasonal  matrices.  The  Ri  are  were  calculated 
according  to  Equation  3.9and  are  not  displayed  either;  only  Pi  and  Gij  for  i  =1-5 
are  displayed. 
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Table  3.1:  Transition  Tables,  T.  Entries  in  tables,  Uj,  are  the  numbers  of  clams  with 
fate  i  at  time  t  +  1  that  were  in  size-class  j  at  time  t.  Time  intervals  were  three 
months.  Spring  was  March  to  June  1995,  summer  was  June  to  September,  fall  was 
September  to  December,  and  winter  was  December  1995  to  March  1996.  Since  the 
fate  of  missing  clams  is  unknown,  n  was  used  instead  of  number  deployed  to  calculate 
parameter  values. _ 
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Table  3.2:  Seasonal  Projection  Matrix  Parameter  Values.  Parameter  values  were 
calculated  by  Equation  3.2.3  from  the  stage-transition  data  in  Table  3.2.3.  Only 
parameters  in  the  lower  right  corner  of  the  8  x  8  matrix  shown  in  Table  3.2  are 
displayed  here:  Pi  (on  the  diagonal)  and  Gij  (on  the  subdiagonal)  for  stages  1-5. 
Since  Gi^  for  stages  A,  B,  and  C  are  all  =  1  and  Ri  vary,  they  are  not  shown  in  these 
matrices, _ 
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Because  of  the  extra  information  afforded  by  using  seasonal  data  on  growth  and 
survival,  we  chose  to  use  these  data,  despite  the  problem  using  them  created  in  in¬ 
corporating  the  reproductive  contribution  data.  The  problem  is  this:  settlement 
measurement  included  survival  to  age  one,  so  the  parameter  measures  reproductive 
contributions  for  an  annual  time  step  rather  than  a  three-month  one.  If  these  param¬ 
eters  were  used  to  quantify  transitions  from  stages  2  through  5  directly  to  stage  1, 
this  would  apply  an  annual  mortality  rate  to  these  offspring  during  their  first  three 
months  of  life,  and  they  would  be  moved  to  stage  1  at  an  age  of  three  months,  rather 
than  one  year.  This  discrepancy  in  timing  was  resolved  by  adding  dummy  stages  A, 
B,  and  C  to  the  model,  as  shown  in  Figure  3.1.  In  our  model,  offspring  are  con¬ 
tributed  to  stage  A  according  to  Ri^  which  include  a  year’s  worth  of  mortality.  But 
then  the  juveniles  pass  through  stages  B  and  C  to  stage  one,  making  each  transition 
with  probability  one.  Therefore,  they  arrive  in  stage  one  at  the  appropriate  survival 
rate  and  at  the  proper  time. 

3.2.4  Estimation  of  Reproduction  Parameters 

To  estimate  the  reproductive  parameters  in  a  matrix  model,  we  need  to  know  size- 
specific  fecundity  and  the  survival  of  those  offspring  to  the  next  time  step.  For  animals 
with  planktonic  larvae,  pre-settlement  survival  is  very  difficult  to  estimate,  and  even 
estimating  post-settlement  survival  challenging,  especially  over  short  time  intervals. 
We  allocated  the  production  of  observed  settling  juveniles  to  size  classes  in  propor¬ 
tion  to  the  reproductive  output  of  that  size  class  (called  “anonymous  reproduction” 
in  Caswell,  1989).  These  estimates  were  adjusted  to  an  annual  interval  by  rough 
estimates  of  survival  rates  from  settlement  to  age  one.  This  method  assumes  that 
the  population  is  closed,  that  is,  that  larvae  spawned  by  these  adults  settle  back  into 
the  same  population.  The  reproductive  contributions  {Ri)  were  calculated  as:  1)  the 
proportion  Mj  of  population  reproductive  output  contributed  by  size  class  i\  2)  the 
number  s  of  recruited  juveniles  per  adult  (by  area);  and  3)  the  probability  S  that 
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settled  juveniles  will  survive  to  one  year  of  age: 

Ri  =  MisS.  (3.9) 

To  estimate  the  Mj,  eggs  per  ripe  female  were  counted  from  histological  prepara¬ 
tions  of  gonad  samples,  using  standard  stereological  methods  (Weibel,  1979).  Game- 
togenesis  begins  in  the  winter  in  Cape  Cod  populations  of  M.  arenaria,  and  spawning 
occurs  throughout  the  summer  from  March  to  September.  Samples  were  taken  from 
five  clams  per  size  class  in  March,  June,  August,  and  September  1995.  No  ripe  fe¬ 
males  in  size  class  two  were  present  in  the  samples  taken.  A  value  was  calculated 
using  numbers  of  eggs  from  clams  collected  in  Quincy,  Massachusetts,  which  were 
nearly  identicle  to  those  of  Barnstable  Harbor  clams  in  other  size  classes.  The  Mj 
(Table  3.2.4)  were  calculated  from  the  mean  number,  Ei,  of  eggs  per  ripe  individual 
in  size  class  i\ 

Reproduction  was  incorporated  into  the  model  in  the  summer  matrix  only.  This 
approach  averaged  reproductive  effort  over  the  spawning  season. 

An  estimate  of  variability  in  settlement  was  obtained  from  data  in  Goshima  (1982). 
He  provides  numbers  of  just-settled  recruits  and  adults  per  area,  at  five  sites  (less 
than  400  m  apart)  over  six  years.  We  estimated  the  ratio  of  newly  settled  clams 
to  adult  clams  (to  nearest  order  of  magnitude)  from  plots  of  clam  density  for  each 
time  and  site  (Goshima,  1982:  Figure  7).  If  the  plotted  value  for  either  young  clam 
density  or  adult  density  was  zero,  we  assumed  that  in  twice  the  area  one  clam  would 
be  found,  and  so  used  0.5  for  our  calculations.  This  was  done  to  avoid  division  by, 
or  taking  the  log  of,  zero.  The  resulting  values  (Table  3.4)  range  over  several  orders 
of  magnitude.  These  numbers  were  divided  by  two  to  count  only  female  offspring  in 
the  model. 

The  survival  probability,  5,  of  the  settlers,  until  the  next  time  step  of  the  model 
(age  one  year),  was  roughly  estimated  as  0.092.  This  value  is  the  mean  of  five  years 
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Table  3.3:  Mean  numbers  of  eggs  per  female  in  size  class  i,  Ei,  sample  size,  n,  the 
variance  of  Ei,  and  proportion  of  population  reproduction  from  size  class  i,  Mi.  Values 
in  parentheses  are  from  Quincy,  Massachusetts,  population  of  clams. 


Size  Class 

n 

Ei 

variance 

Mi 

1 

1 

1.7x10® 

0 

0.02 

2 

(1) 

(3.7x10®) 

0 

0.05 

3 

4 

6.5x10® 

1.2x10^^ 

0.08 

4 

4 

2.03x10^ 

8.7x10^^ 

0.25 

5 

3 

4.77x10^ 

1.8x10^® 

0.60 

Table  3.4:  Numbers  of  newly-settled  clams  per  adult  used  to  characterize  a  distribu¬ 
tion  for  the  stochastic  model,  estimated  from  data  in  Goshima,  1982. 


Site 

1974 

1975 

Year 

1976 

1977 

1978 

1979 

A 

8,000 

0.04 

0.02 

4 

4 

4 

C 

80,000 

0.4 

0.02 

40 

4 

40 

G 

400 

0.4 

0.4 

4 

4 

4 

H 

400 

0.4 

0.002 

4000 

400 

400 

E 

40 

0.04 

0.04 

0.4 

4 

40 
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of  data  on  changes  in  M.  arenaria  cohort  density  with  time,  and  was  estimated  from 
Figure  9  in  Goshima  (1982). 

3.2.5  Population  Structure 

Simulations  were  started  with  the  observed  stage  distribution  in  Barnstable  Harbor. 
This  distribution,  n(0),  was  determined  from  counting  and  measuring  the  clams  col¬ 
lected  in  ten  cores  30  cm  in  diameter  and  30  cm  deep  taken  near  the  study  site  on 
September  26,  1995.  Clams  were  sorted  from  the  sediment  by  hand.  The  population 
consisted  of  3%  class  1,  12%  class  2,  47%  class  3,  28%  class  4,  and  9%  class  5  clams. 

3.3  Model  Analysis  and  Results 

3.3.1  Deterministic  Model 

In  the  deterministic  model,  the  value  used  for  s  was  the  geometric  mean  of  settlement 
data  (Table  3.4),  2.66.  The  deterministic  population  growth  rate,  the  dominant 
eigenvalue  of  A,  was  calculated  as  In  A,  and  equalled  -1.90.  This  population  is 
expected  to  decrease  (Figure  3.2),  but  as  s  is  increased,  population  growth  rate  will 
also  increase(Figure  3.4). 

The  equilibrium  settlement  rate  (ESR)  was  calculated  numerically  by  finding  a 
value  of  s  that  gave  a  In  A  =  0.00.  (An  analytical  method  for  finding  ESR  for  age- 
classified  matrix  models  is  found  in  Vaughan  and  Saila,  1976.)  The  ESR  was  211, 
two  orders  of  magnitude  higher  than  the  geometric  mean  of  values  measured  for  M. 
arenaria  in  Japan  (Goshima,  1982). 

The  stable  stage  distribution  showed  that,  asymptotically,  the  population  will  be 
dominated  by  clams  in  size-class  1  (33%),  with  the  remainder  split  fairly  equally 
between  classes  3,  4,  and  5  (about  20%  each).  This  was  not  the  structure  of  the 
Barnstable  Harbor  population  at  the  time  that  it  was  sampled,  demonstrating  that 


88 


Time  in  Years 


Figure  3.3:  Projections  of  Deterministic  and  Stochastic  Models.  Total  population  size 
is  shown  on  a  logarithmic  scale  through  time  in  years  for  the  deterministic  model  (- 
-  -)  and  ten  randomly  selected  examples  of  simulations  of  the  stochastic  model  (— ). 
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Figure  3.4:  Change  in  deterministic  population  growth  rate,  A,  with  change  in  settle¬ 
ment  rate,  s.  The  *  indicates  the  value  used  in  the  deterministic  model  to  calculate 
the  reported  growth  rate  and  the  o  indicates  the  value  for  the  ESR,  where  the  popu¬ 
lation  is  at  equilibrium.  The  slope  of  this  curve  is  elasticity  of  A  to  s. 
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it  was  not  at  equilibrium. 

Seasonal  elasticities  were  calculated  according  to  the  modification  of  Equation  3.6 
for  periodic  matrices  found  in  Caswell  and  Trevisan  (1994).  Elasticities  show  (Ta¬ 
ble  3.3.1)  that  in  all  seasons,  small  changes  in  survival  in  the  largest  size  class  will  lead 
to  the  largest  change  in  A  of  all  the  parameters.  Elasticity  to  reproduction  increases 
with  size,  but  is  low  overall.  The  slope  of  the  log-log  plot  of  changes  in  A  with  change 
in  s  (Figure  3.3.1)  is  the  elasticity  to  s.  It  increases  as  s  increases. 

Sensitivities  for  the  annual  matix  were  calculated  as: 


1  ^A 
A  dttij 


(3.11) 


for  comparison  to  the  stochastic  sensitivities.  Reproduction  in  the  first  size  class 
had  the  highest  sensitivity  (Table  3.6),  and  reproduction  in  the  larger  size  classes 
and  growth  into  the  largest  size  class  had  high  sensitivities  as  well.  Sensitivities  to 
transitions  from  size  class  two  were  negligible. 


3.3.2  Stochastic  Model 

The  matrix  A  was  varied  stochastically  by  changing  s  at  each  time  step.  All  other 
parameters  were  held  constant.  Values  for  s  for  each  simulation  were  generated 
randomly  (using  Matlab),  from  a  lognormal  distribution.  The  lognormal  distribution 
was  characterized  from  a  normal  distribution  with  mean  and  variance  of  the  logarithm 
of  the  measured  settlement  values  from  Goshima  (1982).  The  anti-logarithm  of  one 
of  these  values  was  then  substituted  into  the  projection  matrix  at  each  time.  An 
example  of  a  typical  time-series  of  settlement  values  generated  by  this  process  shows 
a  stereotypical  random  pattern  if  it  is  plotted  on  a  log  scale,  while  on  a  linear  scale 
,  the  pattern  appears  as  many  years  with  low  settlement  and  a  few  years  with  high 
settlement  (Figure  3.5).  Values  of  s  that  are  higher  than  the  ESR  occur  47%  of  the 
time.  The  skewness  of  the  lognormal  distribution  caused  this  pattern,  which  is  the 
important  consequence  of  using  it. 
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Table  3.5:  Seasonal  elasticity  matrices  for  My  a  arenaria.  Where  seasonal  matrix 
elements  are  zero,  elasticities  are  also  zero.  Those  elements  have  been  replaced  in  the 
displayed  elasticity  matrices  with  dashes,  for  clarity. _ 


Spring 

- 

- 

- 

- 

- 

- 

- 

- 

_ 

0.25 

“ 

- 

- 

- 

- 

- 

- 

- 

<0.01 

<0.01 

- 

- 

- 

- 

- 

- 

<0.01 

0.01 

0.04 

- 

- 

- 

- 

- 

- 

- 

0.15 

0.08 

- 

- 

- 

- 

- 

- 

- 

0.01 

0.44 

Summer 

- 

- 

- 

0.01 

<0.01 

0.02 

0.07 

0.15 

- 

- 

- 

0.04 

<0.01 

- 

- 

- 

- 

- 

- 

0.21 

<0.01 

- 

- 

- 

- 

- 

- 

- 

- 

0.03 

0.03 

- 

- 

- 

- 

- 

- 

- 

0.13 

0.30 

Fall 

. 

. 

.25 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

_ 

- 

- 

- 

<0.01 

- 

- 

- 

- 

- 

- 

- 

<0.01 

0.02 

- 

- 

_ 

- 

- 

- 

- 

0.02 

0.21 

- 

- 

- 

- 

- 

- 

- 

- 

0.06 

- 

- 

- 

- 

- 

- 

- 

0.43 

Summer 

- 

- 

- 

- 

- 

- 

- 

- 

- 

0.25 

_ 

- 

- 

- 

<0.01 

- 

- 

- 

_ 

- 

- 

- 

<0.01 

0.01 

- 

- 

- 

- 

- 

- 

<0.01 

<0.01 

0.19 

- 

- 

- 

- 

- 

- 

- 

0.04 

0.05 

- 

- 

- 

- 

- 

- 

- 

0.01 

0.43 
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Figure  3.5;  Time  Series  of  Stochastic  Settlement.  An  example  of  a  typical  time  series 
of  randomly  generated  settlement  values  plotted  on  a  lognormal  scale  (upper  panel), 
and  on  a  linear  scale  (lower  panel),  showing  the  type  of  pattern  produced  by  the 
stochastic  process. 
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Table  3.6:  Sensitivity  matrix  for  annual  projection  matrix  in  deterministic  model. 


0 

0 

0 

7.08 

0.16 

4.26 

4.85 

4.36 

0 

0 

0 

6.33 

0.15 

3.80 

4.33 

3.89 

0 

0 

0 

2.15 

0.05 

1.29 

1.48 

1.32 

0 

0 

0 

1.70 

0.04 

1.02 

1.17 

1.04 

0 

0 

0 

1.12 

0.03 

0.68 

0.77 

0.69 

0 

0 

0 

0.57 

0.05 

0.34 

0.39 

0.35 

0 

0 

0 

2.28 

0.05 

1.37 

1.57 

1.41 

0 

0 

0 

4.98 

0.12 

2.99 

3.41 

3.07 

Analysis  of  the  stochastic  model  consisted  of  calculating  the  stochastic  growth 
rate,  In  Aj  (Caswell,  1989).  First,  the  one-step  growth  rates  were  calculated  for  each 
time  step  of  a  single  population  trajectory  of  length  T: 

\n  Xs{i)  =  In  N  {i  +  1)  -  \n  N  (i)  i  —  (3.12) 

where  N  is  total  population  size.  Then  the  mean  of  the  one-step  growth  rates  was 
calculated: 

InA,  =  ^:r^^lnA,(i).  (3.13) 

i 

Finally,  the  stochastic  growth  rate  was  reported  as  the  mean  of  the  In  A^  from  many 
individual  trajectories.  In  this  study,  stochastic  growth  rate  was  estimated  as  the 
mean  of  ten  trajectories,  each  run  for  50,000  time  steps  with  the  first  10,000  discarded 
as  transients.  For  the  theoretical  population  of  M.  arenaria  experiencing  variable 
settlement  on  a  lognormal  scale.  In  As  =  —0.95.  Example  of  trajectories  for  the 
stochastic  population  (Figure  3.3.1)  shows  that  it  is  projected  to  decline  less  rapidly 
than  the  deterministic  population. 
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Sensitivities  for  the  stochastic  model, 


dln\ 

ddij 


S 


(3.14) 


were  calculated  according  to  the  stochastic  equivalent  of  Equation  6: 


1  ^ 
s  =  iE 


v,u;_ 


T  V!U,A(i) 


(3.15) 


where  prime  denotes  the  transpose  of  a  matrix,  and  Uf  and  V*  are  the  stochas¬ 
tic  analogs  of  the  deterministic  stable  stage  distribution,  w,  and  reproductive  value 
vector,  V.  These  are  defined  as: 


..  A-t+iUt 

II  A,«U,  II 

(3.16) 

and 

II  • 

(3.17) 

The  symbols  || 

•  II  denote  the  sum  of  the  vector  entries. 

and  the  matrices  At  are 

the  matrices  in  the  sequence  generated  in  the  stochastic  model  (Tuljapurkar,  1990). 
Stochastic  sensitivities  (Table  3.7)  show  the  same  pattern  as  deterministic  ones,  with 
negligible  sensitivity  to  transitions  from  stage  two,  the  highest  sensitivity  to  tran¬ 
sitions  from  stage  one,  and  highest  sensitivity  in  each  column  for  reproduction  and 
growing  into  the  largest  size  classes.  The  stochastic  sensitivities  in  the  upper  left 
of  the  matrix  are  about  double  the  deterministic  ones,  while  in  the  lower  right  they 
are  almost  the  same  as  the  deterministic  ones.  Stochastic  sensitivities  for  transitions 
from  size  class  one  to  size  classes  four  and  five  were  more  than  double  those  in  the  de¬ 
terministic  case,  indicating  that  rapid  growth  to  large  size  would  be  more  important 
to  a  population  with  random  settlement  than  to  one  with  steady  settlement. 

We  were  also  interested  in  knowing  how  strongly  our  results  depended  on  the  esti¬ 
mate  of  variability  in  settlement  that  we  used.  As  a  simple  measure  of  the  sensitivity 
of  stochastic  growth  rate  to  change  in  the  parameters  characterizing  the  lognormal 
settlement  distribution,  we  calculated  the  stochastic  growth  rate  for  ranges  of  these 
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Table  3.7:  Sensitivity  matrix  for  annual  matrix  in  stochastic  model. 


0 

0 

0 

12.73 

0.08 

2.09 

2.55 

1.92 

0 

0 

0 

11.58 

0.08 

1.92 

2.33 

1.76 

0 

0 

0 

4.13 

0.03 

0.69 

0.85 

0.65 

0 

0 

0 

3.50 

0.02 

0.62 

0.76 

0.59 

0 

0 

0 

2.53 

0.02 

0.50 

0.61 

0.48 

0 

0 

0 

1.61 

0.01 

0.38 

0.46 

0.38 

0 

0 

0 

6.01 

0.05 

1.36 

1.66 

1.34 

0 

0 

0 

13.39 

0.12 

3.10 

3.78 

3.06 

values.  We  tested  values  of  the  mean  of  the  log  of  the  settlement  distribution  (holding 
the  variance  constant)  ranging  from  -1.5  to  3,  and  values  of  the  variance  of  the  log 
of  the  settlement  distribution  (holding  the  mean  constant)  ranging  from  0  to  8,  both 
in  increments  of  0.5.  As  the  mean  and  variance  of  the  settlement  distribution  were 
increased,  the  stochastic  growth  rate  increased,  but  at  different  rates  (Figure  3.6). 
For  population  growth  rates  to  reach  equilibrium  (In  =  0),  the  variance  must  be 
almost  doubled  over  the  value  used  at  the  *,  while  the  mean  must  be  more  than 
tripled.  When  the  variance  of  the  distribution  is  equal  to  zero,  the  stochastic  growth 
rate  equals  the  deterministic  growth  rate,  as  expected. 

3.4  Discussion 

Incorporating  variability,  on  a  lognormal  scale,  in  annual  settlement  density  into  a 
matrix  model  for  M.  arenaria  leads  to  a  dramatic  increase  in  population  growth  rate 
compared  to  that  using  mean  settlement  levels.  This  result  suggests  that  variability 
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Mean  Variance 


Figure  3.6;  Sensitivity  of  Stochastic  Growth  Rate  to  Mean  and  Variance  of  Settlement 
Distribution.  The  change  in  stochastic  growth  rate  with  change  in  the  parameters 
specifying  the  distribution  used  in  generating  s  values.  Means  and  variances  were 
calculated  from  the  natural  logs  of  the  settlement  values  in  Table  3.4.  The  change 
in  stochastic  growth  rate  with  change  in  the  mean  of  the  log  of  the  distribution  (left 
panel),  holding  the  variance  constant  at  the  value  (4.20)  used  to  calculate  the  reported 
A.  The  change  in  stochastic  growth  rate  with  change  in  the  variance  of  the  log  of 
the  distribution,  holding  the  mean  constant  at  the  value  (0.98)  used  to  calculate  the 
reported  A.  The  *  indicates  the  value  of  the  mean  (left  panel)  or  variance  (right 
panel)  used  to  calculate  the  reported  A. 


97 


in  larval  settlement  is  a  critical  factor  determining  population  state.  Observed  large 
standing  stocks  of  M.  arenaria  may  be  due  to  one  good  year  of  set,  but  the  popula¬ 
tion  persists  because  clams  live  long  enough  that  another  good  set  is  likley  to  occur 
before  they  all  die.  Previous  analyses  of  deterministic  models  (Brousseau  et  ai,  1982; 
Malinowski  and  Whitlatch,  1988)  have  likely  underestimated  the  potential  for  growth 
in  M.  arenaria  populations. 

These  results  contrast  with  those  of  Gotelli  (1991),  who  concluded  that  the  con¬ 
tribution  of  variability  in  recruitment  to  temporal  dynamics  of  a  broadcast-spawning 
coral  was  minimal.  He  created  the  variability  in  his  model  by  randomizing  twenty- 
three  monthly  matrices  in  series.  Variability  in  recruitment  (measured  as  new  colonies 
appearing  in  the  first  size  class,  each  month)  ranged  over  only  three  orders  of  magni¬ 
tude  (from  0.00  to  1.22;  pers.  comm.,  N.  Gotelli),  as  opposed  to  the  eight  orders  of 
magnitude  of  variability  found  in  our  model.  This  might  be  due  to  a  shorter  plank¬ 
tonic  period  (3-20  d)  for  the  coral  than  for  M.  arenaria  (14-21  d).  The  increased 
amount  of  variability  in  our  model  could  easily  explain  why  we  found  it  to  be  im¬ 
portant  but  Gotelli  did  not,  especially  since  we  also  found  that  increased  variability 
increased  stochastic  growth  rate  (Figure  3.6). 

Since  the  model  is  sensitive  to  the  mean  and  variance  of  the  distribution  used 
(Figure  3.6),  uncertainty  as  to  the  parameters  describing  this  distribution  would  ob¬ 
viously  lead  to  uncertainty  in  the  population  growth  rate.  Since  we  do  not  know  if 
the  mean  and  variance  of  the  settlement  distribution  used  here  are  similar  to  actual 
settlement  in  Barnstable  Harbor,  there  is  no  reason  to  conclude  that  the  population 
there  is  expected  to  collapse.  In  fact,  because  our  results  show  that  the  theoretical 
clam  population  is  decreasing  so  rapidly,  we  believe  there  is  some  mismatch  between 
growth  and  survival  in  Barnstable  Harbor  and  the  settlement  levels  in  Japan.  The 
ESR  calculated  was  two  orders  of  magnitude  higher  than  mean  settlement  in  Japan, 
but  that  population  was  not  in  decline  (Goshima,  1982),  so  clams  must  have  higher 
survival  rates  there  than  in  Barnstable,  MA.  Settlement  in  Barnstable  Harbor  mea- 
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sured  in  September,  1995  (Ripley  and  McDowell,  in  prep.)  was  0.43.  This  means  that 
about  63%  of  settlement  events  in  our  model  were  larger  than  measured  settlement 
in  Barnstable,  suggesting  that  settlement  ranges  in  Japan  may  be  similar  to  those  in 
Barnstable  Harbor. 

Equilibrium  settlement  rates  have  been  used  to  estimate  larval  survival,  under  the 
assumption  that  the  population  is  in  equilibrium  (Vaughan  and  Saila,  1976;  Brousseau 
et  ai,  1982).  In  this  study,  the  population  was  not  in  equilibrium  when  it  was  surveyed 
(the  population  structure  was  not  the  same  as  the  stable  stage  distribution),  and  we 
believe  that  few  field  populations  of  marine  bivalves  are  in  equilibrium,  making  this 
assumption  highly  quesitonable.  The  ESR  calculated  was  much  higher  than  observed 
settlement,  yet  the  field  population  of  clams  in  Barnstable  have  been  a  productive 
resource  for  decades.  Given  the  importance  of  variation  in  settlement  demonstrated 
here,  we  believe  that  calculation  of  ESR  is  of  little  informational  value. 

Elasticities  calculated  here  for  reproductive  parameters  in  the  deterministic  case 
were  low,  which  agree  with  Gotelli’s  results.  Further,  Gotelli  (1991)  found  in  a  random 
simulation  of  1000  stage-structured  matrices  that  elasticities  for  survivorship  were 
always  larger  than  those  for  growth,  and  that  reproductive  elasticities  were  always 
smallest.  However,  since  they  are  defined  over  “small”  changes  in  parameter  values, 
elasticities  tell  us  little  about  the  importance  of  large  changes  in  the  parameters. 
The  results  of  this  study  showed  that  despite  small  deterministic  elasticities,  random 
change  in  reproductive  parameters  led  to  fairly  large  change  in  population  growth 
rate.  Nakaoka  (1997)  also  found  that  the  recruitment  parameters  had  low  elasticities 
but  could  cause  large  changes  in  A  when  they  were  made  stochastic. 

Variation  in  growth  rates  or  survival,  certainly  plays  a  role  in  population  pro¬ 
cesses  of  M.  arenaria  which  we  did  not  address  in  this  study.  However,  variability  in 
recruitment  was  shown  to  play  a  larger  role  in  determining  population  growth  than 
variation  in  shell  growth  rates  by  Nakaoka  (1997).  Indeed,  since  growth  and  survival 
parameters  are  probabilities  and  must  be  between  zero  and  one,  it  is  always  possible 
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for  the  fecundity  parameters  to  vary  much  more  widely.  Growth  rates  are  also  more 
likely  to  vary  over  much  smaller  ranges,  which  would  cause  less  dramatic  changes 
in  population  growth  rates.  It  is  important  to  select  the  appropriate  vital  rates  to 
vary  in  a  model  to  answer  speciific  questions  about  their  role  in  population  processes 
(Noda  and  Nakao,  1996;  Nakaoka,  1997),  and  here  we  chose  to  focus  on  settlement. 

The  settlement  data  used  here  are  lognormally  distributed  (according  to  a  normal 
probability  plot;  not  shown),  which  means  that  the  distribution  is  skewed,  with  few 
large  values  and  many  more  small  ones.  Other  marine  invertebrate  recruitment  data 
sets  have  also  shown  lognormal  trends.  The  distribution  Nakaoka  (1997)  used  was 
lognormal,  estimated  from  a  large  data  set  (ca.  20  years)  of  recruitment  in  Yoldia 
notabilis,  a  bivalve  which  is  similar  in  life  history  to  M.  arenaria.  Noda  and  Nakao 
(1996)  observed  massive  recruitment  for  a  species  of  marine  gastropod  once  in  ev¬ 
ery  eight  years,  with  little  or  none  in  between,  also  suggesting  a  lognormal  pattern. 
Recruitment  patterns  for  marine  invertebrates  have  repeatedly  been  shown  to  have 
this  kind  of  pattern  {e.g.  Coe,  1956;  Loosanoff,  1964;  Caffey,  1985),  so  our  results 
may  be  broadly  applicable.  Further  research  is  needed  in  characterizing  recruitment 
dynamics  of  intertidal  organisms  so  that  it  is  possible  to  describe  year-to-year  and 
spatial  variability  more  completely.  It  would  also  be  interesting  to  investigate  density- 
dependence  of  recruitment,  which  might  damp  the  changes  in  population  size  due  to 
changes  in  settlement  generated  in  this  density-independent  model.  Characterizing 
how  recruitment  varies  with  different  developmental  modes  will  also  be  crucial  to 
understanding  marine  invertebrate  population  dynamics. 

What  is  the  role  of  variability  in  settlement  in  the  life  history  of  marine  inverte¬ 
brates?  Results  here  suggest  that  the  advantage  of  broadcasting  larvae  may  not  be 
solely  in  greater  dispersal  (Palmer  and  Strathmann,  1981;  Strathmann  and  Strath- 
mann,  1982),  but  also  in  periodic  large  recruitment  successes.  For  animals  with  a  life 
span  long  enough  to  sustain  populations  through  several  years  with  no  settlement, 
broadcast  spawning  is  a  gamble  with  potentially  big  payoffs.  Alternatively,  the  low 
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parental  investment  of  energy  into  planktonic  larvae,  which  leads  to  variability  in 
settlement,  may  have  led  to  the  possibility  of  long  life  span.  In  either  case,  adult  life 
span  and  variable  recruitment  are  linked. 

3.5  Summary 

1.  Equilibrium  settlement  rates  provide  little,  if  any,  information  on  settlement  impact 
on  population  structure. 

2.  Variability  in  settlement  is  apparently  a  gamble  worth  taking  for  My  a  arenaria 
populations.  This  conclusion  is  expected  to  hold  true  for  other  organisms  with  high 
fecundity,  uncertain  recruitment,  and  relatively  long  life  spans. 

3.  While  deterministic  models  may  be  useful  in  comparing  populations  to  one  an¬ 
other,  stochastic  models  present  a  much  different  picture  of  population  processes  for 
organisms  with  highly  variable  vital  rates,  if  that  variability  can  be  characterized. 
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Chapter  4 


The  interaction  of  chronic 
contaminant  exposure  and 
ecological  factors  in  population 
processes  of  an  infaunal  clam 

4.1  Introduction 

The  physiological  effects  of  lipophilic  organic  contaminants  on  marine  bivalve  molluscs 
have  been  examined  extensively  during  the  past  decade.  The  majority  of  the  studies 
have  been  conducted  on  the  blue  mussel  Mytilus  edulis  {e.g.,  Bayne  et  al,  1985; 
Bayne  et  al,  1988)  with  an  effort  to  integrate  responses  over  several  levels  of  biological 
hierarchy  and  to  examine  responses  linked  to  specific  classes  of  contaminants.  The 
general  results  of  this  research  are  that  tissue  burdens  of  lipophilic  organic  compounds 
may  interfere  with  normal  metabolic  processes  influencing  growth,  development,  and 
reproduction  (Capuzzo  et  al.,  1988).  Exposure  to  contaminants  is  known  to  impair 
physiological  mechanisms  (Capuzzo  and  Sasner,  1977;  Gilfillan  et  al,  1977,  Widdows, 
1985),  co-occur  with  histopathological  disorders  (Moore,  1988;  Lowe,  1988),  and 
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cause  loss  of  reproductive  potential  (Berthou  et  ai,  1987;  Neff  and  Haensly,  1982). 
Recent  work  has  extended  this  approach  to  other  species  of  bivalve  molluscs,  such  as 
the  subtropical  turkey  wing  mussel,  Area  zebra  (Addison  and  Clarke,  1990;  Widdows 
et  ai,  1990)  and  the  soft  shell  clam,  Mya  arenaria  (Leavitt  et  ai,  1990;  Weinberg  et 
al.,  1997;  McDowell  Capuzzo  et  ai,  in  prep.). 

The  soft-shell  clam,  Mya  arenaria,  is  common  and  commercially  important  in 
temperate  mudflats  along  the  eastern  US.  It  lives  at  high  densities  at  both  clean  and 
contaminated  sites,  such  as  New  Bedford  Harbor,  MA,  and  the  Boston  Harbor,  MA, 
area.  Even  in  areas  where  Mytilus  edulis  has  oocyte  atresion  to  the  point  that  they 
are  not  reproducing,  M.  arenaria  appears  to  be  reproducing  normally  (McDowell 
Capuzzo  and  Leavitt,  1995).  Gardner  and  Pruell  (1988),  however,  found  significant 
histopathological  lesions  in  soft  shell  clams  at  selected  contaminated  sites  in  Quincy 
Bay  including  gill  inflammation,  atypical  cell  hyperplasia  in  gill  and  kidney,  hyper¬ 
parasitism  with  rickettsia  in  digestive  ducts/tubules,  and  general  parasitism.  Similar 
results  were  found  by  Moore  et  al.  (1996). 

Changes  in  population  structure  and  dynamics  are  likely  results  of  alterations  to 
reproductive  and  developmental  potential,  which  are  expected  results  of  chronic  expo¬ 
sure  to  chemical  contaminants.  Alterations  in  bioenergetics  linked  with  observations 
of  reduced  fecundity  and  viability  of  larvae,  abnormalities  in  gamete  and  embryolog- 
ical  development,  and  reduced  reproductive  success  provide  a  strong  empirical  basis 
for  examination  of  population  responses.  Koojiman  and  Metz  (1984)  suggested  that 
such  sublethal  effects  of  contaminant  exposure  should  be  interpreted  in  light  of  the 
survival  probabilities  and  reproductive  success  of  populations,  bridging  the  gap  be¬ 
tween  individual  and  population  responses.  Although  a  wide  range  of  sublethal  stress 
indices  have  been  proposed  for  evaluation  of  chronic  responses  of  organisms  to  con¬ 
taminants,  few  have  been  linked  to  the  survival  potential  of  the  individual  organism 
or  the  reproductive  potential  of  the  population. 

From  the  known  deleterious  effects  of  contaminant  exposure  on  bivalves,  we  would 
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expect  M.  arenaria  living  at  contaminated  sites  to  be  unhealthy.  However,  we  find 
dense  populations  of  this  clam  at  contaminated  sites.  The  observed  population  toler¬ 
ance  to  chemical  stress  may  be  due  to  one  or  more  of  the  following  reasons;  1)  healthy 
larvae  spawned  at  clean  sites  recruiting  to  contaminated  sites  and  subsidizing  popula¬ 
tions;  2)  reductions  in  survival,  growth,  or  reproduction  are  not  sufficient  to  suppress 
population  growth;  3)  predators  (including  humans)  are  less  tolerant  than  clams  to 
chemical  contaminant  exposure;  4)  organic  enrichment  at  contaminated  sites  provides 
nutrition  for  clams;  and  5)  adaptation  to  chemical  contaminant  exposure.  As  a  first 
step  in  understanding  the  differences  in  populations  found  at  contaminated  and  un¬ 
contaminated  sites,  we  assessed  (2)  by  measuring  clam  growth,  survival,  fecunditity, 
and  recruitment  at  five  sites  across  a  gradient  in  contaminant  concentrations.  Using  a 
matrix  population  model,  vital  rates  were  integrated  to  evaluate  the  population  level 
effects  of  differences  in  individual  physiological  parameters. 

The  matrix  population  model  can  be  thought  of  as  an  analytical  tool  that  com¬ 
bines  vital  rates  of  individuals  into  information  about  the  population  as  a  whole.  If 
we  assume  that  potential  biochemical,  cellular,  and  physiological  changes  associated 
with  contaminant  stress  are  manifested  as  changes  at  the  organismal  level,  namely  in 
fecundity,  survival,  and  growth,  we  can  say  that  the  matrix  model  integrates  these 
factors  and  addresses  their  combined  effect  at  the  population  level. 

To  construct  a  matrix  population  model,  the  population  is  divided  into  categories 
that  correspond  to  changes  in  vital  rates  over  the  life  span  of  the  animal.  In  this 
case,  size-classes  were  used,  since  for  M.  arenaria,  important  life  history  traits  such 
as  predation  rates  (Blundon  and  Kennedy,  1982)  and  fecundity  (Brousseau,  1978)  are 
size-specific,  and  because  size  is  more  easily  measured  than  age.  The  set  of  linear 
equations  that  describe  the  proportion  of  clams  that  move  between  these  categories 
in  a  given  time  step,  are  written,  in  matrix  notation,  as: 

n(t-f- 1)  =  An(t).  (4.1) 

The  vectors  n  are  size-class  frequency  distributions  and  A  is  the  projection  matrix. 
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Each  element,  (the  matrix  element  in  the  ith  row  and  jth  column  of  A),  represents 
contributions  from  class  j  to  class  i  over  one  time  step.  Elements  on  the  diagonal 
of  A  represent  survival  and  stasis,  and  are  noted  Pi.  Elements  on  the  subdiagonal 
represent  survival  and  growth  between  size-classes,  and  are  noted  Gij.  Contributions 
to  the  smallest  size  class  via  reproduction  appear  in  the  first  row  of  the  matrix,  and 
are  noted  here  as  jR,. 

Analysis  of  the  model  includes  calculation  of  the  dominant  eigenvalue  of  A,  A, 
which  is  the  asymptotic  population  growth  rate  and  a  measure  of  mean  fitness  of 
individuals  in  the  population  (Caswell,  1989).  Elasticities  are  the  proportional  sensi¬ 
tivities  of  A  to  each  matrix  element,  and  are  the  relative  contributions  of  the  matrix 
elements  to  fitness.  These  contributions  can  be  compared  between  populations  by 
doing  a  life  table  response  experiment  (LTRE)  decomposition  (Caswell,  1996).  By 
treating  the  whole  life  table  as  the  response  variable,  one  can  calculate  where  in  the 
life  table  the  differences  in  A  between  populations  are  found  (Caswell,  1989).  Elastic¬ 
ity  analysis  and  LTRE  decompositions  are  powerful  methods  to  evaluate  what  parts 
of  the  life  cycle  have  the  most  impact  on  population  processes  {e.g.  Levin  et  ai,  1996, 
and  other  work  cited  in  Caswell,  1996). 

4.2  Methods 

4.2.1  Study  Sites 

Based  on  surveys  of  lipophilic  organic  contaminant  concentrations  and  histopathological 
conditions  of  soft  shell  clam  populations,  we  selected  five  sites  from  those  previously 
studied  in  Boston  Harbor  and  Massachusetts  Bay  (Moore  et  al,  1996;  Shea  and 
Seavey,  1994).  Three  contaminated  sites  around  Boston — Saugus  River,  Eort  Point 
Channel,  and  Neponset  River — and  two  reference  sites  about  100  km  away  from 
Boston  on  Cape  Cod — Wellfleet  and  Barnstable  Harbors — were  selected  (Figure  4.1). 
On  the  basis  of  sedimentary  PAH  concentrations,  these  sites  reflect  a  gradient  of 
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contamination  and  observed  histopathological  effects  in  soft  shell  clam  populations. 
Sites  were  similar  in  sediment  type  (sandy  mud). 

Sediment  and  clam  tissue  contaminant  concentrations  (Table  4.1)  were  evaluated 
during  the  sampling  for  this  study  (McDowell  and  Shea,  1996).  Clams  for  chemical 
analyses  were  collected  with  the  clams  used  for  analyses  described  in  this  study,  and 
sediment  samples  for  chemical  analysis  were  taken  at  the  locations  where  clams  were 
deployed  in  the  mark-recapture  study.  Low  levels  of  contaminants  are  present  even 
at  Barnstable  and  Wellfleet  Harbors,  the  Neponset  site  is  moderately  contaminated, 
and  the  Saugus  River  and  Fort  Point  Channel  sites  are  highly  contaminated.  The 
differences  between  sediment  and  clam  concentrations  of  these  chemicals  reflects  dif¬ 
ferences  in  bioavailablity  and  bioaccumulation  of  specific  compounds,  and  perhaps 
saturation  threshholds,  as  discussed  in  McDowell  and  Shea  (1996). 

4.2.2  Adult  Growth  and  Survival 

Growth  and  survival  probabilities  were  measured  in  a  mark-recapture  study  (similar 
to  Weinberg  et  ai,  1996)  at  each  site.  Measurements  were  made  separately  during 
four  periods  of  time  over  the  year,  as  growth  and  survival  of  M.  arenaria  are  known 
to  change  seasonally  (Brousseau,  1978,  1979).  In  March  1995,  sufficient  numbers  of 
clams  were  dug  from  the  intertidal  mudflat  to  select  twenty  in  each  of  the  following 
shell  length  size  classes:  20-40  mm  (class  1) ;  40-49.9  mm  (class  2)  ;  50-59.9  mm  (class 
3)  ;  60-69.9  mm  (class  4)  ;  and  >70  mm  (class  5).  Size  classes  were  selected  to  include 
all  sizes  present  in  the  populations.  Although  clam  growth  slows  with  size,  clams  were 
expected  to  grow  enough  to  measure  transitions  between  larger  size-classes. 

Each  clam  was  measured  and  marked  with  an  identification  number  on  both  valves 
in  indelible  ink.  At  this  time,  clams  at  each  site  were  numbered  from  one  to  one 
hundred.  Ten  clams  (two  from  each  size  class)  were  placed  into  separate  holes  poked 
in  each  of  ten  sediment-filled  mesh  bags  (mesh  size  ca.  20  mm)  buried  so  that  the 
surface  of  the  sediment  that  they  contained  was  flush  with  the  surrounding  sediment. 
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Figure  4.1:  Map  of  study  sites  in  Massachusetts,  USA.  Latitudes  are  north,  longitudes 


Table  4.1;  Sediment  and  Clam  Tissue  Contaminant  Concentrations  (ng/g  dry  weight). 
Polycyclic  aromatic  hydrocarbons  (PAH),  polychlorinated  biphenyls  (PCB),  dioxins 
(DDT),  and  linear  alkyl  benzenes  (LAB)  in  sediment  samples  and  in  clam  tissues. 
Values  are  the  means  of  two  samples  each  of  sediments  and  clams  collected  in  March, 
1995  (from  McDowell  and  Shea,  1996).  Clams  were  selected  randomly  from  the  pop¬ 
ulation  sampled  for  the  mark-recapture  study.  Sediments  were  collected  no  farther 
than  1  m  from  locations  where  clams  were  deployed  in  the  mark-recapture  study. 


Site 

Total  PAH 

sed.  clams 

Total  PCB 

sed.  clams 

Total  DDT 

sed.  clams 

Total  LAB 

sed.  clams 

Barnstable  H. 

352 

300 

0.5 

11.5 

0.2 

9.1 

7 

80 

Wellfleet  H. 

102 

367 

<0.1 

14.7 

0.2 

0.4 

10 

360 

Neponset  R. 

1,450 

1,900 

5.4 

130.5 

2.0 

14.4 

31 

1,680 

Saugus  R. 

18,342 

5,110 

12.8 

91.1 

5.1 

20.6 

65 

910 

Fort  Point  C. 

66,121 

7,370 

34.8 

56.7 

26.8 

14.4 

925 

8,400 
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in  their  normal  orientation  with  siphon  up  and  at  the  appropriate  depth  for  their  size 
(about  one  body  length).  Bags  were  located  near  the  center  of  the  region  of  the  mud 
flats  exposed  at  mean  low  water.  The  bags  did  not  exclude  predators,  but  simply 
marked  the  location  of  the  clams  and  allowed  for  easier  recovery.  The  density  of  clams 
within  the  bags  (470/m^)was  about  half  the  mean  of  natural  densities  at  the  sites 
(820/m^).  In  June,  September,  and  December,  the  bags  were  removed,  the  clams 
were  sorted  from  the  sediment,  and  growth  and  mortality  rates  determined  by  noting 
the  numbers  and  valve  lengths  of  the  clams  found  alive  or  dead.  At  each  sampling 
time,  new  clams  were  collected  to  replace  missing  or  dead  clams,  labelled,  and  mixed 
with  the  survivors  before  placing  all  clams  in  new  bags  as  in  March  1995.  New 
numbers  were  used  to  label  the  clams  at  each  time,  so  that  they  were  not  confused 
with  survivors  from  the  previous  time  period  (he.  clams  first  deployed  in  June  were 
numbered  starting  with  101;  those  first  deployed  in  September  were  numbered  starting 
with  200;  and  those  first  deployed  in  December  were  numbered  starting  with  300). 
Surviving  clams  that  were  re-deployed  were  assumed  to  have  the  same  survival  and 
growth  rates  as  newly  deployed  clams.  In  March  1996,  all  clams  were  removed  from 
the  bags,  counted,  and  measured. 

An  additional  set  of  ten  bags  was  deployed  at  each  site  in  March  1995  and  not 
removed  until  March  1996,  to  assess  the  effect  of  regularly  disturbing  clams  in  the 
three-month  deployments. 

There  were  clams  missing  from  the  bags  on  most  sampling  dates,  and  on  occasion 
entire  bags  were  missing.  All  of  the  winter  bags  in  Saugus  were  missing  in  March 
1996.  To  obtain  data  on  winter  growth  and  survival  in  Saugus,  another  set  of  bags 
was  deployed  at  this  site  during  the  next  year.  These  bags  were  in  place  from  January 
to  April  1997.  Some  difference  in  growth  and  survival  between  years  is  expected,  so 
these  data  were  not  the  best  to  use  to  compare  to  the  other  sites  in  a  different  year, 
but  they  are  the  only  data  available. 

Missing  clams  were  probably  either  washed  out  of  bags  after  they  had  died,  or 
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were  removed  by  predators.  We  found  some  of  the  valves  near  the  bags  on  the  surface 
of  the  sediment,  so  some  of  the  missing  clams  were  clearly  dead.  They  could  not  all 
be  dead,  however,  because  when  we  assumed  that  they  were  dead,  the  model  results 
predicted  that  the  probability  of  finding  large  clams  would  be  much  smaller  than  the 
actual  probability  of  finding  them  in  the  field.  Due  to  this  uncertainty  over  the  fate 
of  missing  clams,  we  based  calculation  of  matrix  parameters  only  on  the  clams  that 
were  found.  This  assumes  that  the  proportion  of  live  clams  is  the  same  for  the  missing 
clams  as  for  the  ones  we  found.  If  this  assumption  is  wrong,  parameter  estimates  are 
either  low  or  high,  depending  on  whether  a  smaller  or  larger  proportion  of  the  missing 
clams  are  alive  than  was  the  case  for  the  clams  we  found. 

Matrix  elements  Pi  and  Gij  for  i,j  =  1-5  were  calculated  for  each  size  class, 
season,  and  site  from  the  field  data.  The  stasis  parameters.  Pi,  and  were  calculated 
as  the  number  of  clams  in  size  class  i  at'  season  t  that  were  found  alive  and  in  size 
class  i  at  season  t  +  1,  divided  by  the  total  number  of  clams  in  size  class  i  at  season 
t  that  were  found  (alive  or  dead)  at  season  t  +  1.  The  growth  parameters,  Gij,  were 
calculated  as  the  number  of  clams  at  size  class  i  at  season  t  that  were  found  alive  in 
size  class  j  a,t  t  +  1,  divided  by  the  total  number  of  clams  in  size  class  i  at  season  t 
that  were  found  (alive  or  dead)  at  season  t  +  1. 

4.2.3  Reproductive  Cycle  and  Fecundity 

To  document  the  annual  reproductive  cycle  and  measure  fecundity,  five  clams  from 
each  of  the  five  size  classes  selected  for  growth  and  survival  measurements  were  col¬ 
lected  from  the  intertidal  area  at  each  site  in  March,  June,  August,  Septmeber,  and 
December.  The  clams  were  brought  into  the  laboratory  where  they  were  measured, 
dissected,  and  the  digestive  gland-gonad  (DG)  complex  was  weighed.  Aliquots  of 
gonad  tissue  were  cut  from  the  DG  complex,  weighed,  and  fixed  in  10%  formalin  in 
0.45  fim  filtered  seawater.  Samples  were  embedded  in  paraffin,  sliced  into  ca.  8  fxm 
sections,  mounted  on  slides,  and  stained  with  hematoxylin-eosin  (Humason,  1972). 
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Prepared  slides  from  each  clam  were  examined  for  gender  and  reproductive  con¬ 
dition  according  to  the  definitions  of  indifferent,  early  developing,  late  developing, 
ripe,  spawning,  and  spent  gonadal  condition  in  Mya  arenaria  (Coe  and  Turner,  1938; 
Ropes  and  Stickney,  1965).  To  estimate  fecundity,  eggs  were  counted  from  prepared 
slides  for  female  clams  in  a  late  developing  or  ripe  stage,  using  standard  stereological 
techniques  (Weibel,  1979)  that  have  been  used  previously  to  estimate  fecundity  of  M. 
arenaria  {e.g.  Brousseau,  1978).  Estimates  were  made  from  examination  of  slides 
under  a  light  microscope,  using  a  digitizing  pad.  Sigma  Scan  software,  and  a  camera 
lucida  to  count  and  measure  eggs.  The  mean  number  of  nucleated  eggs  per  unit  area 
was  calculated  for  one  gonad  section  per  clam,  based  on  counts  of  nucleated  eggs  in 
10  non-overlapping  10x10  unit  retucule  grids.  The  mean  oocyte  diameter  and  nuclear 
diameter  of  25  nucleated  eggs  per  gonad  section  were  measured.  The  number  of  eggs 
per  unit  volume  of  gonad  were  calculated  by  converting  grid  units  to  area  (mm^),  di¬ 
viding  by  the  number  of  grid  squares  (100),  and  dividing  by  mean  nuclear  diameter. 
The  total  number  of  eggs  per  gonad  was  calculated  based  on  the  total  volume  of  the 
DG  complex  (calculated  from  data  on  weight-to- volume  ratio  of  the  DG  complex),  the 
total  percent  of  gonadal  tissue  within  the  digestive  gland-gonad  complex  (D.  Leavitt, 
unpublished  data),  and  the  number  of  eggs  per  volume  of  gonadal  tissue. 

4.2.4  Recruitment  and  Population  Structure 

Reproductive  contributions  in  matrix  models  (the  first  row  of  A)  are  usually  denoted 
Fi  and  are  estimated  as  fecundity  multiplied  by  the  probablity  that  offspring  survive 
to  the  next  model  time  step.  This  method  is  problematic  for  soft  shell  clams  and 
other  marine  invertebrates  that  have  planktonic  larvae,  because  we  know  so  little 
about  transport  patterns  and  survival  during  both  the  planktonic  period  and  early 
post-settlement  life.  One  can  estimate  fecundity  from  histological  sections  of  gonads, 
but  it  is  difficult  to  estimate  larval  survival.  Matrix  models  have  actually  been  used 
to  estimate  larval  survival  (Brousseau  et  al.,  1982;  Vaughan  and  Saila,  1976),  but 
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these  models  cannot  also  be  used  to  obtain  a  population  growth  rate  because  it  is 
necessary  to  assign  a  value  to  A  to  calculate  larval  survival. 

In  this  study,  we  defined  recruits  as  clams  that  had  settled  and  survived  at  the 
study  site  to  age  one  year.  The  ratio  of  recruits  to  adults,  r,  is  used  in  the  model  as  the 
basis  of  the  reproductive  contribution  parameter.  This  method  avoids  the  problem 
of  not  being  able  to  quantify  larval  survival  because  it  matches  the  time  scale  over 
which  we  can  keep  track  of  young  clams.  The  parameter  used  here  is  denoted  Ri  to 
emphasize  this  different  method  of  estimation.  Ri  are  calculated  as; 

Ri  =  O.brMi  (4.2) 

where  Mi  is  the  proportion  of  population  reproductive  output  contributed  by  size 
class  i.  This  was  calculated  as  the  mean  fecundity  of  females  in  class  i,  divided 
by  the  sum  of  mean  fecundities  over  all  size  classes.  Using  this  factor  maintains 
the  population  structure  of  the  model,  but  assumes  that  the  population  is  closed 
with  respect  to  recruitment.  This  method  also  assumes  that  offspring  from  different 
sizes  of  clams  have  equal  survival  probabilities,  and  that  the  population  that  actually 
spawned  the  settling  juveniles  (not  necessarily  the  one  at  the  study  site)  has  the  same 
size  structure  as  the  study  population.  Recruitment  is  multiplied  by  0.5  to  count  only 
females,  assuming  a  1:1  sex  ratio  (Brousseau,  1978). 

To  measure  r,  we  had  to  be  able  to  identify  one-year  (365  day)  old  clams.  We 
estimated  the  size  of  one  year  old  clams  by  measuring  more  than  200  clams,  collected 
in  August  1995,  from  Barnstable  Harbor,  known  to  have  set  on  a  flat  with  no  other 
clams  on  it  about  one  year  prior  to  sampling  (pers.  com.,  Tom  Marcotti,  Barnstable 
Shellfish  Biologist).  These  clams  ranged  in  size  from  17  to  40  mm,  with  a  mean  of 
26.49  ±  3.98.  The  upper  limit  of  this  range  fell  so  closely  to  the  size  class  cut-off 
of  40  mm  that  clams  less  than  that  size  were  assumed  to  be  one  year  old  or  less. 
This  assumption  doubtless  underestimates  recruitment  in  Barnstable  Harbor,  since 
clams  that  settled  in  the  spring  will  be  larger  than  this  size,  but  may  overestimate 
recruitment  at  the  other  sites  if  those  clams  grow  more  slowly. 
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Recruitment  and  clam  densitites  in  the  field  were  estimated  at  each  site  in  ten 
replicate  30  cm-diameter  by  30  cm-deep  cores  by  measuring  and  counting  the  clams. 
Cores  were  taken  haphazardly  in  areas  where  clams  were  present,  as  evidenced  by 
their  siphon  holes  at  the  sediment  surface.  Cores  were  taken  in  September  1995, 
about  one  year  from  the  last  spawning  event  of  the  previous  year,  and  again  in  June 
1996.  Clams  were  picked  from  the  sediment  by  hand,  measured,  and  counted.  The 
ratio  of  clams  less  than  40  mm  to  those  larger  than  40  mm  (r)  was  calculated  for  each 
site  based  on  total  numbers  measured  in  1995. 

4.2.5  Data  Analysis 

The  sites  we  studied  differ  in  many  factors  besides  contaminant  concentrations,  so 
they  cannot  be  treated  as  representing  individual  contaminant  levels.  Hypotheses 
about  differences  between  clams  at  uncontaminated  and  contaminated  sites  were 
tested  by  treating  the  two  least  contaminated  sites,  Barnstable  and  Wellfieet,  as 
replicates  of  the  clean  condition,  and  the  two  most  contaminated  sites,  Fort  Point 
Channel  and  Saugus,  as  replicates  of  the  “contaminated”  condition.  Everywhere  that 
these  comparisons  are  made,  the  words  “clean”  and  “contaminated”  appear  in  quota¬ 
tion  marks  to  alert  the  reader.  Hypotheses  concerning  recruitment,  growth,  fecundity, 
and  survival  were  tested  in  this  way.  Note  that  these  are  not  interpreted  as  tests  for 
the  effects  of  contaminants  on  clams,  merely  as  tests  for  differences  between  clams 
living  at  sites  where  there  are  contaminants  present  and  those  where  contaminants 
are  not  present. 

Where  statistical  tests  were  done  using  data  based  on  size  classes,  the  original  size 
classes  were  used,  not  the  revised  classes  (!’  to  5’)  used  to  calculate  matrix  paramters. 
All  means  are  reported  i  one  standard  error.  ANOVAs  and  t— tests  were  performed 
using  the  SPSS  statistical  package  (v.  6.1,  SPSS  Inc.,  Chicago,  IL),  and  log-linear 
anlyses  were  performed  using  Systat  (v.  7.0,  SPSS  Inc.  Chicago,  IL). 
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4.2.6  Model  Formulation  and  Analysis 

Because  of  the  extra  information  afforded  by  using  seasonal  data  on  growth  and 
survival,  we  chose  to  use  these  data,  despite  the  disparity  in  time  step  length  between 
these  data  and  reproductive  contribution  data.  The  problem  is  this:  recruitment  was 
measured  as  settlement  and  survival  to  age  one,  so  the  parameter  is  over  an  annual 
time  step  rather  than  a  three-month  step,  as  in  the  case  of  growth  and  survival  data. 
If  these  parameters  were  used  to  quantify  transitions  from  stages  2  through  5  directly 
to  stage  1,  this  would  apply  an  annual  mortality  rate  to  these  offspring  during  their 
first  three  months  of  life,  and  they  would  be  moved  to  stage  1  at  an  age  of  three 
months,  rather  than  one  year.  This  discrepancy  in  timing  was  resolved  by  adding 
dummy  stages  A,  B,  and  C  to  the  model  (Figure  4.2).  In  our  model,  offspring  are 
contributed  to  stage  A  according  to  Ri,  which  include  a  year’s  worth  of  mortality.  But 
then  the  juveniles  pass  through  stages  B  and  C  to  stage  one,  making  each  transition 
with  probability  one.  Therefore,  they  arrive  in  stage  one  at  the  appropriate  survival 
rate  and  at  the  proper  time. 

Parameters  for  seasonal  matrices  at  each  site  were  calculated  as  described  above. 
Seasonal  matrices,  A^,  were  multiplied  together  in  sequence  {A-^inter  x  A  fail  x  Asummer 
X  Aspring)  to  yield  an  annual  projection  matrix,  A,  for  each  site.  As  it  turned  out, 
larger  clams  were  growing  too  slowly  to  grow  between  size  classes  at  all  sites  but  Barn¬ 
stable  Harbor.  This  caused  the  annual  matrices  at  these  sites  to  be  reducible,  that  is, 
it  was  not  possible  for  clams  to  move  into  some  stages.  Irreducibility  is  required  for  all 
of  the  properties  of  eigenvalues  and  eignevectors  described  by  the  Perron-Frobenius 
theorem,  and  the  subsequent  analysis  of  the  model,  to  hold  (Caswell,  1989). 

To  remedy  this  problem,  the  data  were  regrouped  into  a  set  of  size  classes  which 
were  narrower  in  width  as  the  clams  got  larger:  20-40  mm  (class  1’);  40-54.9  mm 
(class  2’);  55-59.9  mm  (class  3’);  60-64.9  mm  (class  4’);  and  >65  mm  (class  5’).  This 
change  caused  there  to  be  fewer  clams  starting  in  size  classes  3’  and  4’  than  in  3  and 
4.  Although  the  new  size  classes  did  not  make  the  Wellfleet,  Saugus,  Neponset,  or 
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Figure  4.2;  Life  Cycle  Diagram  for  Mya  arenaria.  The  circles  indicate  stages,  and  the 
arrows  represent  these  transitions:  surviving  and  growing  to  a  subsequent  stage  (Gy), 
surviving  and  staying  in  the  same  stage  (Fj),  or  reproductive  contribution  to  the  first 
stage  (Fj).  Stages  A,  B,  and  C  are  dummy  stages,  designed  to  align  the  timing  of 
recruitment  parameters  with  growth  and  survival  parameters.  The  corresponding 
seasonal  projection  matrix  is  found  in  Table  4.2. 
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Table  4.2:  Elements  of  the  seasonal  projection  matrices,  A*,  corresponding  to  the 
stage  transitions  in  Figure  4.2.6.  All  locations  in  the  matrix  where  no  parameter 
is  shown  are  zero.  In  this  model,  G^b,c,  and  Gc,i  are  always  one  and  are 

not  displayed  in  subsequent  tables  of  seasonal  matrices.  Since  the  Ri  only  occurr  in 
the  spring  or  summer  matrices  (depending  on  site)  and  are  calculated  according  to 
Equation  4.2,  they  are  not  displayed  in  subsequent  seasonal  matrix  tables  either;  only 
Pi  and  Gij  for  z=l-5  are  displayed. 
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Fort  Point  Channel  matrices  reducibile,  it  reduced  the  number  of  changes  required. 
The  Wellfleet  Harbor  annual  matrix  was  made  reducible  by  replacing  the  zero- valued 
6*1,2  parameter  in  the  summer  matrix  with  0.01  (an  arbitrarily  chosen  small  value). 
Saugus  and  Neponset  annual  matrices  were  made  reducible  by  adding  0.01  to  the 
zero-valued  63,4  parameter  in  the  summer  matrix.  Only  two  clams  were  found  alive 
at  Fort  Point  Channel  in  the  fall,  so  zero-valued  parameters  P3,  P4,  61,2,  63,4,  and 
64,5  had  to  be  replaced  with  0.01  to  make  the  annual  matrix  reducible. 

Population  growth  rates,  and  elasticities  for  periodic  matrices  were  calculated 
according  to  Caswell  (1989)  and  Caswell  and  Trevisan  (1994). 

Since  recruitment  rates  are  known  to  be  variable  (Goshima,  1982),  some  investiga¬ 
tion  of  how  large  changes  (rather  than  the  small  ones  implicit  in  analytical  sensitivity 
analysis)  in  the  recruitment  parameter  might  effect  the  population  growth  rate  was 
desirable.  This  was  tested  numerically  by  calculating  population  growth  rates  using 
the  highest  and  lowest  vaues  of  recruitment  measured  for  any  site,  to  compare  with 
population  growth  rates  calculated  with  the  measured  site-specific  recruitment  rates. 

We  calculated  an  LTRE  decomposition  matrix  to  determine  which  matrix  param¬ 
eters  caused  the  difference  between  A  values  at  the  “clean”  and  “contaminated”  sites. 
The  difference  in  population  growth  rate,  A  A,  between  the  annual  matrix  cal¬ 
culated  from  the  means  of  the  seasonal  matrices  from  Fort  Point  Channel  and  Saugus 
River,  and  the  matrix  the  annual  matrix  calculated  from  the  means  of  the 

seasonal  matrices  from  Barnstable  and  Wellfleet  Harbors  was  approximated  by: 


EE(^5&) 


where  M 


+  A^ 


where  s  are  seasons,  bij  are  elements  of  the  seasonal  matrices  for  “clean”  and  “con¬ 
taminated”  sites,  and  A  are  annual  matrices  for  these  sites.  Essentially,  the  peri¬ 
odic  sensitivities  of  A  to  the  mean  of  clean  and  contaminated  matrices  is  multiplied 
element-by-element  by  the  element-by-element  difference  between  the  two  matrices 
(Caswell,  1989),  for  each  season. 
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4.3  Results 


4.3.1  Adult  Growth  and  Survival 

It  is  likely  that  survival  and  growth  of  clams  in  the  field  experiment  were  compromised 
by  the  stress  caused  by  digging  them  up,  drying  and  handling  them,  and  replacing 
them  in  the  sediment.  An  assessment  of  the  impact  of  handling  on  clam  growth  can 
be  obtained  from  comparing  growth  rates  of  clams  in  the  one-year  study  to  the  clams 
in  the  three-month  studies.  For  clams  in  the  three-month  bags,  growth  over  the  year 
was  calculated  as  difference  in  length  between  March  1995  and  1996.  The  majority 
of  clams  surviving  the  year  in  the  one-year  bags  came  from  Saugus  and  Neponset 
(69/71),  while  in  the  three-month  bags,  the  majority  (11/15)  came  from  Neponset. 
As  no  meaningful  comparison  could  be  made  between  sites,  all  sites  were  pooled. 
Differences  in  growth  between  inter-sampling  periods  was  tested  with  size-class  in  a 
2-way  ANOVA  (Table  4.3).  No  difference  was  detected  between  growth  in  one  year 
(Af  =  4.44  ±  0.75  mm)  and  the  four  three— month  [X  =  4.13  ±  2.0  mm)  treatments 
p=0.67). 

More  clams  were  found  at  the  end  of  the  one-year  treatments  than  at  the  end  of 
the  fourth  three-month  ones.  More  clams  may  have  washed  out  of  the  three-month 
bags  since  the  sediments  in  the  bags  left  undisturbed  for  a  year  were  more  compact 
and  stable  than  those  which  had  been  reworked  during  sampling  three  times.  Growth 
rates  did  not  differ  between  the  lengths  of  time  bags  were  deployed,  and  since  clams 
were  treated  the  same  at  each  site,  we  feel  that  handling  stress  is  not  an  important 
factor  to  consider  in  making  comparisons  between  these  sites. 

Growth  rates  were  only  calculated  for  clams  found  alive.  Differences  in  clam 
growth  between  “clean”  and  “contaminated”  sites  was  tested  against  size  class  and 
season  by  3— way  ANOVA  (Table  4.4).  All  three  factors  were  highly  significant,  as 
well  as  all  interactions  except  site— by— size  class.  Mean  growth  over  all  seasons  and 
size  classes  was  1.9  ±0.15  mm/3  months  for  clams  at  the  “clean”  site  and  0.42  ±0.11 
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Table  4.3;  ANOVA  table  for  growth  in  one-year  and  three-month  bags.  Source  of 
variation  (Source)  is  explained  by  the  sum  of  squares  (SS),  degrees  of  freedom  (df), 
mean  square  (MS),  F  statistic  (F),  and  significance  (p).  “Timing”  is  the  variable  for 
length  of  time  bags  were  deployed. _ 


Source 

SS 

df 

MS 

F 

P 

Main  Effects 

1206.1 

5 

241.2 

10.7 

<0.001 

Timing 

4.2 

1 

4.2 

0.2 

0.67 

Size  Class 

1205.1 

4 

301.4 

13.3 

<0.001 

Timing  x  Size  Class 

31.0 

4 

7.8 

0.3 

0.85 

Explained 

1882.4 

9 

209.2 

9.3 

<0.001 

Residual 

1717.8 

76 

22.6 

Total 

3600.2 

85 

42.4 

mm/3  months  at  the  “contaminated”  site.  Further,  growth  was  higher  at  the  “clean” 
site  any  way  means  were  compared.  Clams  in  larger  size  classes  grew  much  more 
slowly  than  smaller  clams  (Figure  4.3). 

Affects  of  site,  season,  and  size-class  on  survival  were  compared  by  log-linear 
analysis,  since  survival  data  are  frequencies.  Number  of  clams  in  each  size  class  at 
time  t  that  were  dead  or  alive  at  time  t  + 1  were  tallied.  Under  the  null  hypothesis  of 
fate  (being  dead  or  alive)  being  dependent  on  initial  state  and  independent  of  site  and 
season,  models  including  state  and  fate  but  with  different  combinations  of  site,  season, 
and  interaction  terms  were  compared.  Significance  was  tested  on  the  difference  of 
and  degrees  of  freedom  between  models  which  were  identical  except  for  the  term(s) 
in  question  (Caswell,  1989).  Models  are  hierarchical  and  are  specified  here  by  the 
highest  level  interaction  terms  that  they  include.  Model  comparisons  (Table  4.5) 
showed  highly  significant  site  (FT,  SFT)  and  season  (FL,  SFL)  effects.  These  effects 
were  still  highly  significant  when  tested  on  models  that  included  the  other  factor. 
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Table  4.4:  ANOVA  Table  for  growth  rates.  Source  of  variation  (Source)  is  explained 
by  the  sum  of  squares  (SS),  degrees  of  freedom  (df),  mean  square  (MS),  F  statistic 
(F),  and  significance  (p). 


Source 

SS 

df 

MS 

F 

P 

Main  Effects 

3620.9 

8 

452.6 

64.5 

<0.001 

Site 

408.6 

1 

408.6 

58.2 

<0.0001 

Size 

1530.1 

4 

382.5 

54.5 

<0.001 

Season 

1263.9 

3 

421.3 

60.0 

<0.001 

2-Way  Interactions 

865.5 

19 

45.5 

6.5 

<0.001 

Site  X  Size 

46.2 

4 

11.5 

1.6 

0.16 

Site  X  Season 

304.2 

3 

101.4 

14.4 

<0.0001 

Size  X  Season 

504.4 

12 

42.0 

5.9 

<0.0001 

Site  X  Size  x  Season 

169.57 

12 

14.1 

2.0 

0.02 

Explained 

4641.1 

39 

119.0 

16.9 

<0.0001 

Residual 

7018.9 

1000 

7.0 

Total 

11660.3 

1039 

11.2 
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Figure  4.3:  Mean  growth  rates  of  clams  in  size  classes  1-5  for  “clean”  (open  bars)  and 
“contaminated”  (filled  bars)  sites.  Means  are  over  individuals  and  seasons.  Negative 
growth  rates  are  due  to  normal  wear  and  chipping  of  shells  of  larger  clams,  which  do 
not  grow  as  rapidly  as  their  shells  are  worn  down.  Error  bars  represent  one  standard 
error. 
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showing  that  they  each  explain  significant  amounts  of  the  difference  in  fates.  The 
interaction  of  time  and  location  (FTL,  SFTL)  was  also  highly  significant.  The  mean 
(over  sites  and  times)  of  survival  frequencies  was  higher  at  the  “contaminated”  site 
(0.83)  than  at  the  “clean”  site  (0.54),  and  most  of  this  difference  occurred  in  the 
summer  and  fall  (Figure  4.4). 

For  a  general  picture  of  differences  in  growth  and  survival  between  all  the  sites,  we 
compared  multiway  contingency  tables  of  initial  state  and  subsequent  fate  of  clams 
by  general  log-linear  model  fitting.  This  analysis  was  identical  to  that  for  survival, 
except  that  fates  were  not  dead/alive,  but  were:  size  classes  F-fi’,  dead,  or  missing. 
Contingency  table  analysis  results  (Table  4.6)  show  highly  significant  season  and  site 
effects,  even  when  both  factors  are  included  in  the  model.  This  is  evidence  that 
seasonal  matrices  calculated  from  these  contingency  tables  are  significantly  different 
from  one  another. 

Sites  were  also  compared  with  log-linear  analysis  separately  by  season  to  try  to 
break  down  the  effects  of  the  full  multiway  tables.  This  simpler  analysis  tested  the 
saturated  (SFL)  against  the  SL,  SF  model  to  find  the  significance  of  the  FL,  SFL 
model.  There  were  highly  significant  differences  between  sites  in  all  seasons  but 
winter  (Table  4.7). 

Growth  rates  across  the  entire  size  spectrum  at  each  site,  for  each  season,  were 
compared  in  Walford  plots  of  length  at  time  t  against  length  at  time  t  +  1  (Wal- 
ford,  1946).  Sites  were  compared  for  difference  in  slope  of  linear  regressions  through 
Walford  plots  by  analysis  of  covariance  of  length  data  (McCuaig  and  Green,  1983). 
Tests  for  the  homogeneity  of  slope  for  regression  lines  of  Walford  plots  for  each  season 
(Figure  4.5)  were  significant  in  all  seasons,  indicating  that  the  slopes  did  differ  (Table 
4.8).  However,  in  fall  and  winter  (but  not  spring  or  summer),  if  Barnstable  Harbor  is 
removed  from  the  analysis,  the  other  sites  do  not  differ  significantly  from  one  another, 
suggesting  that  the  slope  of  the  Barnstable  Harbor  line  is  different  enough  from  the 
others  to  cause  the  non-significant  result. 
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Table  4.5:  Results  of  Log-Linear  Analysis  of  Effects  of  Site  (L),  Season  (T),  and 
Size  Class  (S)  on  Fate  (F).  Fates  are  being  dead  or  alive.  The  two  models  that  were 
compared  to  calculate  and  degrees  of  freedom  for  each  set  of  factors  being  tested 
are  shown  above  the  horizontal  lines. 


Model 

G2 

df 

P 

STL,  SF 

512.5 

35 

STL,  SFT 

151.8 

20 

FT,  SFT 

360.7 

15 

<0.0001 

STL,  SF 

512.5 

35 

STL,  SFL 

384.9 

30 

FL,  SFL 

127.6 

5 

<0.0001 

STL,  SFT 

151.8 

20 

STL,  SFT,  SFL 

36.5 

15 

SFT,  FT 

115.3 

5 

<0.0001 

STL,  SFL 

384.9 

30 

STL,  SFL,  SFT 

36.5 

15 

SFT,  FT 

348.4 

15 

<0.0001 

STL,  SFL,  SFT 

36.5 

15 

SFTL 

0 

0 

FTL,  SFTL 

36.5 

15 

0.001 
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Figure  4.4:  Survival  frequencies  by  season.  Means  are  over  individuals  in  all  size 
classes.  Open  bars  are  for  the  “clean”  site  and  filled  bars  are  the  “contaminated” 
site,  for  each  season. 
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Table  4.6;  Results  of  Log-Linear  Analysis  of  Effects  of  Site  (L),  Season  (T),  and  Size 
Class  (S)  on  Fate  (F).  Fates  are  the  five  size  classes,  dead,  or  missing.  The  two  models 
that  were  compared  to  calculate  and  degrees  of  freedom  for  each  set  of  factors 
being  tested  are  shown  above  the  horizontal  lines. 


Model 

G2 

df 

P 

STL,  SF 

1342.9 

570 

STL,  SFT 

842.1 

480 

FT,  SFT 

500.8 

90 

<0.0001 

STL,  SF 

1342.9 

570 

STL,  SFL 

917.5 

450 

FL,  SFL 

425.4 

120 

<0.001 

STL,  SFT 

842.1 

480 

STL,  SFT,  SFL 

373.0 

360 

FL,  SFL 

469.1 

120 

<0.0001 

STL,  SFL 

917.5 

450 

STL,  SFL,  SFT 

373.0 

360 

FT,  SFT 

544.5 

90 

<0.0001 

STL,  SFL,  SFT 

373.0 

360 

SFTL 

0 

0 

FTL,  SFTL 

373.0 

360 

0.31 
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Figure  4.5:  Walford  Plots  of  Clam  Growth  by  Season.  Length  (in  mm)  at  deployment 
(Start  Size)  is  plotted  against  length  at  recovery  three  months  later  (End  Size)  for  all 
clams  recovered  alive  in  each  season.  Lines  are  regressions  of  the  form  y  =  a*x  +  b 
for  clams  at  each  site.  Regression  parameters  and  comparisons  between  slopes  are 
reported  in  Table  4.3.1. 
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Table  4.7:  Log-linear  analysis  of  contingency  tables  of  fate  (F)  by  state  (S)  in  each 
season.  Fates  are  the  five  size  classes,  dead,  or  missing.  For  each  season,  the  SL, 
SF  model  was  compared  to  the  saturated  model  to  test  for  location  effects.  The 
table  shows  log-liklihood  ratios,  degrees  of  freedom,  and  significance  for  the  FL,  SFL 
model.  _ 


Season 

df 

P 

Spring 

211.6 

120 

<0.001 

Summer 

331.7 

120 

<0.001 

Fall 

204.3 

120 

<0.001 

Winter 

94.5 

120 

0.96 

Interpretation  of  these  results  is  complicated  by  the  possibility  of  having  large 
enough  samples  (40-80  clams  in  each  regression)  that  no  difference  is  detectected 
when  in  fact  there  is  one  (McCuaig  and  Green,  1983).  Walford  plots  are  normally 
made  from  measurements  of  annual  bands  on  bivalve  shells,  whereas  here,  we  used 
sequential  measurements  over  three-month  intervals.  The  time  interval  should  not 
matter,  but  on  our  plots,  each  point  is  an  independent  measurement  on  an  individual 
clam,  rather  than  several  points  being  from  the  same  clam.  These  results  suggest 
that  at  least  in  the  fall  and  winter,  small  clams  in  Barnstable  Harbor  are  growing 
faster  in  proportion  to  large  clams  than  are  small  clams  at  the  other  sites. 

4.3.2  Reproductive  Cycle  and  Fecundity 

Reproductive  cycles  of  clam  populations  differed  between  contaminated  and  uncon¬ 
taminated  sites  (Figure  4.6).  Both  female  and  male  clams  from  Barnstable  Harbor 
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Table  4.8:  Walford  Plot  Regression  Parameters  for  each  site  in  each  season,  and 
the  significance  of  the  test  for  homogeneity  of  slope.  Parameters  are  the  slope,  a, 
and  intercept,  b,  from  regressions  of  the  form  y  =  a  *  x  +  b.  Significance  values  in 
parentheses  are  for  tests  where  Barnstable  Harbor  was  removed  from  analysis. 


Season 

Site 

Slope 

Intercept 

r2 

Significance 

Spring 

Barnstable  H. 

0.81 

15.87 

0.96 

p<0.001 

Fort  Point  C. 

0.85 

9.39 

0.97 

Neponset  R. 

0.93 

4.61 

0.98 

Saugus  R. 

0.86 

8.76 

0.98 

Wellfleet  H. 

0.94 

6.22 

0.94 

Summer 

Barnstable  H. 

0.88 

10.12 

0.94 

p<0.001 

Fort  Point  C. 

0.74 

17.21 

0.92 

Neponset  R. 

0.77 

15.42 

0.95 

Saugus  R. 

0.90 

6.72 

0.98 

Wellfleet  H. 

0.94 

4.33 

0.97 

Fall 

Barnstable  H. 

0.81 

13.23 

0.92 

p<0.001 

Fort  Point  C. 

0.99 

-0.30 

0.99 

(p=0.49) 

Neponset  R. 

0.96 

2.00 

0.99 

Saugus  R. 

0.95 

2.22 

0.99 

Wellfleet  H. 

0.97 

2.33 

0.99 

Winter 

Barnstable  H. 

0.90 

0.93 

0.96 

p<0.001 

Fort  Point  C. 

0.98 

0.87 

0.99 

(p=0.79) 

Neponset  R. 

0.98 

0.64 

0.99 

Saugus  R. 

0.98 

0.60 

0.99 

Wellfleet  H. 

0.98 

0.50 

0.99 

132 


and  Wellfleet  showed  evidence  of  advanced  stages  of  gamete  development  and  spawn¬ 
ing  during  the  late  spring  through  early  fall.  Populations  from  the  Boston  area  (Fort 
Point  Channel,  Saugus  and  Neponset  River)  did  not  show  evidence  of  spawning  until 
mid-summer  and  spawning  occurred  for  only  a  short  period  of  time. 

Only  51  ripe  females  from  all  sites  and  times  were  collected  in  this  sampling 
program.  Because  data  points  were  missing  for  one  or  more  size  classes  at  each  site, 
size-specific  fecundities  (Mj)  could  not  be  calculated  for  each  site.  Pooled  data  from 
all  sites  and  seasons  was  used  to  estimate  Mi  for  each  size  class  1-5  as  0.03,  0.08,  0.12, 
0.27,  and  0.52.  Fecundities  were  not  normally  distributed,  so  egg  numbers  were  log- 
transformed  to  produce  normal  distributions  before  analysis.  Differences  in  fecundity 
between  “clean”  and  “contaminated”  sites  by  size  class  and  sampling  month  were 
tested  by  3-way  ANOVA.  Sampling  month  was  not  significant  (p=0.22),  and  cases 
were  missing,  so  it  was  removed  as  a  factor,  and  site  was  tested  by  size  class  in  a 
2-way  ANOVA  (Table  4.9).  No  difference  was  detected  between  eggs  per  female  of 
“clean”  {X  =  15.98  ±0.19,  n=31)  and  “contaminated”  {X  =  15.71  ±0.19,  n=14) 
sites  (p=0.37).  The  site  by  size  interaction  was  not  significant  (p=0.63). 

4.3.3  Recruitment  and  Population  Structure 

Populations  of  clams  at  the  five  sites  showed  differing  patterns  of  size-class  frequencies 
and  population  density  in  1995  and  1996  (Figure  4.7).  The  Neponset  River  site  had 
the  highest  density  in  both  years.  Recriutment  by  year  at  “clean”  and  “contaminated” 
sites  was  tested  by  2-way  ANOVA.  Recruits  (clams<  40  mm)  per  area  did  not  differ 
significantly  at  the  0.05  level  between  “clean”  {X  =  1.8  ±  .57)  and  “contaminated” 
(X  =  0.95  ±  0.15)  sites(Table  4.10).  Recruits  per  area  did  differ  between  years 
(p=0.01),  and  the  interaction  between  site  and  year  was  also  significant  (p<0.001). 
Recruitment  by  year  was  also  tested  by  2-way  ANOVA  as  recruits  per  adult,  since 
this  was  calculated  for  model  parameter  estimates.  Recruits  per  adult  (r  in  the 
model)  differed  significantly  between  “clean”  {X  =  0.46  ±0.18)  and  “contaminated” 
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Figure  4.6:  Reproductive  Cycle  of  Mya  arenaria  at  each  site  over  a  one-year  period. 
Months  of  the  year  are  labelled  from  1  (January)  to  12  (December).  Number  of  clams 
out  of  samples  (n=25)  collected  in  March,  June,  August,  September,  and  December 
1995  in  a  reproductive  condition  (late  developing,  .ripe,  or  spawning)  are  shown  for 
each  site.  Open  symbols  are  uncontaminated  sites  and  filled  symbols  are  contaminated 
sites.  Sites  are  Barnstable  Harbor  (open  circles).  Fort  Point  Channel  (filled  triangles), 
Neponset  River  (filled  circles),  Saugus  River  (filled  squares),  and  Wellfleet  Harbor 
(open  squares). 
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Table  4.9:  ANOVA  table  for  fecundity. Source  of  variation  (Source)  is  explained  by 
the  sum  of  squares  (SS),  degrees  of  freedom  (df),  mean  square  (MS),  F  statistic  (F), 
and  significance  (p). 


Source  of  Variation 

SS 

df 

MS 

F 

Sig.  ofF 

Main  Effects 

17.2 

5 

3.4 

10.2 

<0.001 

Site 

0.3 

1 

0.3 

0.8 

0.37 

Size 

17.1 

4 

4.3 

12.8 

<0.001 

Site  X  Size 

0.8 

4 

0.2 

0.6 

0.63 

Explained 

27.4 

9 

3.0 

9.1 

<0.001 

Residual 

11.7 

35 

0.3 

Total 

39.1 

44 

0.9 
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—  0.15  ±  0.03)  sites  at  the  0.001  level  (Table  4.10).  Year  and  the  site-by-year 
interaction  were  also  significant  (p<0.001).  There  are  fewer  degrees  of  freedom  for  this 
test  because  recruits  per  adult  could  not  be  calculated  in  samples  where  there  were 
no  adults.  Values  for  recruits  per  adult  depend  on  the  density  of  adult  populations, 
and  since  adult  densities  are  higher  at  contaminated  sites  but  raw  recruitment  does 
not  differ,  recruitment  per  adult  is  higher  at  clean  sites. 

4.3.4  Model  Analysis  Results 

According  to  results  of  the  histological  examination  of  gonad  samples,  clams  in  Barn¬ 
stable  and  Wellfleet  Harbors  spawned  throughout  the  spring  and  summer,  whereas 
clams  from  the  Boston  area  sites  reproduced  during  a  short  period  of  time  in  late 
summer.  This  pattern  was  approximated  in  the  model  by  having  reproduction  occur 
during  spring  and  summer  for  Barnstable  and  Wellfleet,  and  only  during  the  summer 
at  the  Boston  sites.  At  these  times,  the  elements  of  the  first  row  of  the  seasonal  ma¬ 
trix  were  calculated  by  Equation  4.2,  using  measured  values  of  Mi  (Section  4.3.2)  and 
r  (Table  4.18).  =  0  in  seasons  when  reproduction  is  not  occurring.  Only  the  lower 

right  portion  (Table  4.2)  of  seasonal  matrices  containing  Pi  and  Gij  for  i,j  =1-5  are 
reported  (Table  4.11,  Table  4.12,  Table  4.13,  Table  4.14,  and  Table  4.15),  since  Gij 
for  stages  A,  B,  and  C  are  always  one. 

At  each  site,  elasticities  to  each  parameter  were  consistent  across  the  seasons.  The 
largest  element  of  each  matrix  (Table  4.16)  contributed  from  74-100%  to  A.  All  of 
these  parameters  were  stasis  parameters  (Pj),  for  the  largest  size  class  at  clean  sites, 
and  for  the  second  or  third  size  class  at  the  other  sites. 

The  annual  projection  matrices  (Table  4.3.4)  summarize  the  differences  among  the 
sites  more  clearly  than  do  the  seasonal  matrices.  Parameters  in  these  matrices  are 
calculated  from  all  the  ways  clams  can  make  each  transitions  during  a  year.  Matrices 
with  more  entries  below  the  diagonal  reflect  faster  growth  between  size  classes,  and 
the  sums  of  columns  reflect  survival  rates  by  size  class.  Clam  populations  at  Wellfleet 
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Figure  4.7:  Clam  Population  Structures  in  1995  (upper  panel)  and  1996  (lower  panel). 
For  each  site  on  the  x-axis — BH  (Barnstable  Harbor),  FP  (Fort  Point  Channel),  QN 
(Neponset  River),  SR  (Saugus  River),  and  WH  (Wellfleet  Harbor) — the  five  bars  show 
the  mean  number  of  clams  in  ten  30  cm  diameter  cores  for  each  size  class:  1  (open 
bars);  2  (vertical  lines);  3  (diagonal  lines);  4  (horizontal  lines),  and  5  (filled  bars). 
Error  bars  represent  one  standard  error  over  the  ten  replicate  cores. 


137 


Table  4.10:  ANOVA  tables  for  testing  raw  recruitment  (clams<40  mm/core  (30 
diameter))  and  recruits  per  adult. 


Source  of  Variation 

SS 

df 

MS 

F 

Sig.  of  F 

Raw  Recruitment 

Main  Effects 

51.4 

2 

25.7 

4.6 

0.013 

Year 

37.8 

1 

37.8 

6.8 

0.011 

Site 

13.6 

1 

13.6 

2.4 

0.122 

Year  x  Site 

74.1 

1 

74.1 

13.3 

<0.001 

Explained 

125.5 

3 

41.8 

7.5 

<0.001 

Residual 

423.4 

76 

5.6 

Total 

548.9 

79 

6.9 

Recruits  per  Adult 

Main  Effects 

11.9 

2 

5.9 

33.9 

<0.001 

Year 

8.2 

1 

8.2 

46.4 

<0.001 

Site 

6.7 

1 

6.7 

37.9 

<0.001 

Year  x  Site 

9.6 

1 

9.6 

54.5 

<0.001 

Explained 

15.8 

3 

5.2 

29.9 

<0.001 

Residual 

10.2 

58 

0.18 

Total 

26.1 

61 

0.4 

cm 
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Table  4.11:  Barnstable  Harbor  seasonal  projection  matrices 


Winter 

0.71 

0 

0 

0 

0 

0.286 

0.92 

0 

0 

0 

0 

0 

0.818 

0 

0 

0 

0 

0.182 

1 

0 

0 

0 

0 

0 

0.933 

Spring 

0 

0 

0 

0 

0 

1 

0.333 

0 

0 

0 

0 

0.407 

0.143 

0 

0 

0 

0.074 

0.429 

0.25 

0 

0 

0 

0.429 

0.25 

0.867 

Summer 

0 

0 

0 

0 

0 

0.1 

0.143 

0 

0 

0 

0 

0.0476 

0 

0 

0 

0 

0.0476 

0.1 

0 

0 

0 

0 

0.1 

0 

0.333 

Fall 

0.111 

0 

0 

0 

0 

0.333 

0.286 

0 

0 

0 

0 

0.0952 

0.4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.333 
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Table  4.12:  Fort  Point  Channel  seasonal  projection  matrices.  Arbitrary  parameters 
added  to  make  annual  matrix  reducible  are  in  bold. 


Winter 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0.75 

0 

0 

0 

0 

0 

0.8 

Spring 

0.667 

0 

0 

0 

0 

0.333 

0.733 

0 

0 

0 

0 

0.133 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0.786 

Summer 

0.5 

0 

0 

0 

0 

0.25 

1 

0 

0 

0 

0 

0 

0.5 

0 

0 

0 

0 

0.5 

0.6 

0 

0 

0 

0 

0 

0.5 

Fall 

0 

0 

0 

0 

0 

0.01 

1 

0 

0 

0 

0 

0 

0.01 

0 

0 

0 

0 

0.01 

0.01 

0 

0 

0 

0 

0.01 

0.0526 
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Table  4.13:  Neponset  River  seasonal  projection  matrices.  Arbitrary  parameters  added 
to  make  annual  matrix  reducible  are  in  bold. 


Winter 

0.833 

0 

0 

0 

0 

0.0833 

0.944 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0.75 

0 

0 

0 

0 

0.75 

0 

0 

0 

0 

0.25 

0.733 

Spring 

0.917 

0 

0 

0 

0 

0.083 

0.958 

0 

0 

0 

0 

0 

0.727 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0.556 

Summer 

0.222 

0 

0 

0 

0 

0.778 

0.773 

0 

0 

0 

0 

0.091 

0.75 

0 

0 

0 

0 

0.01 

0.5 

0 

0 

0 

0 

0 

0.682 

Fall 

0.667 

0 

0 

0 

0 

0.25 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0.636 
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Table  4.14:  Saugus  River  seasonal  projection  matrices.  Arbitrary  parameters  added 
to  make  annual  matrix  reducible  are  in  bold.  Parameters  for  winter  matrix  were 
estimated  from  data  collected  in  1996.  All  other  data  was  collected  in  1995. 


Winter  1  0 

0  1 

0  0 

0  0 

0  0 

Spring  0.688  0 

0.313  0.885 

0  0.0769 

0  0 

0  0 

Summer  0.167  0 

0.667  0.937 

0  0.0625 

0  0 

0  0 

Fall  0.833  0 

0.167  1 

0  0 

0  0 

0  0 


0  0  0 

0  0  0 

1  0  0 

0  0.83  0 

0  0.083  0.96 

0  0  0 

0  0  0 

1  0  0 

0  1  0 

0  0  0.93 

0  0  0 

0  0  0 

1  0  0 

0.01  0.75  0 

0  0.25  0.619 

0  0  0 

0  0  0 

1  0  0 

0  1  0 

0  0  0.958 
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Table  4.15:  Wellfleet  Harbor  seasonal  projection  matrices.  Arbitrary  parameters 
added  to  make  annual  matrix  reducible  are  in  bold. 


Winter 

1 

0 

0 

0 

0 

0 

0.963 

0 

0 

0 

0 

0 

0.917 

0 

0 

0 

0 

0.083 

0.889 

0 

0 

0 

0 

0 

0.972 

Spring 

0.077 

0 

0 

0 

0 

0.308 

0.259 

0 

0 

0 

0 

0.185 

0.5 

0 

0 

0 

0.074 

0.5 

0.6 

0 

0 

0 

0 

0.4 

0.75 

Summer 

0 

0 

0 

0 

0 

0.01 

0 

0 

0 

0 

0 

0.034 

0 

0 

0 

0 

0 

0.2 

0 

0 

0 

0 

0 

0.083 

0.087 

Fall 

0 

0 

0 

0 

0 

0.091 

0.055 

0 

0 

0 

0 

0 

0.4 

0 

0 

0 

0 

0.2 

0.2 

0 

0 

0 

0 

0 

0.217 
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Table  4.16:  Seasonal  Elasticities.  Only  the  largest  element  of  elasticity  matrices  for 
each  season  and  site  are  reported  since  the  other  elements  are  negligible  in  compari¬ 
son.  Barnstable  Harbor(BH),  Fort  Point  Channel(FP),  Neponset  River  (QN),  Saugus 
River  (SR),  and  Wellfleet  Harbor  (WH). 


Site 


BH 

FP 

QN 

SR 

WH 

Element 

Ps 

P2 

P2 

Pz 

Pz 

Spring 

0.89 

0.99 

0.75 

0.97 

1 

Summer 

0.90 

0.99 

0.83 

0.97 

1 

Fall 

0.91 

0.99 

0.74 

0.97 

1 

Winter 

0.91 

0.99 

0.75 

0.97 

1 
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are  clearly  limited  by  survival,  while  growth  rates  are  highest  at  Barnstable  Harbor. 
Survival  rates  are  highest  at  Neponset  River  and  Saugus.  When  A  is  calculated  from 
the  annual  matrices,  Barnstable  and  Wellfleet  Harbors  have  the  lowest  population 
growth  rates  (Table  4.3.4).  The  highest  elasticities  in  seasonal  matrices  were  to 
stasis  parameters.  Pi,  which  suggests  that  survival  may  be  important  for  M.  arnearia 
populations.  The  relative  magnitude  of  A  is  robust  to  changes  in  r  over  several  orders 
of  magnitude,  and  in  fact  does  not  change  population  growth  rate  by  more  than 
0.01  even  when  varied  over  a  wide  range  (Table  4.3.4),  as  is  expected  given  the  low 
elasticity  to  these  parameters.  Population  growth  is  higher  at  contaminated  sites 
than  at  uncontaminated  ones  (Figure  4.8). 

Population  growth  rate  for  the  “clean”  site  matrix  was  0.046,  and  for  the  “con¬ 
taminated”  site  matrix  it  was  0.79,  which  are  close  to  the  means  of  the  population 
growth  rates  of  the  matrices  used  to  calculate  the  “clean”  (Barnstable  and  Wellfleet) 
and  “contaminated”  (Fort  Point  and  Saugus)  site  matrices.  The  largest  elements 
of  the  seasonal  LTRE  decomposition  matrices,  Dj,  (Table  4.19)  are  those  which  are 
principally  responsible  for  the  observed  difference  in  population  growth  rates  between 
the  two  populations.  By  far  the  largest  element  in  each  seasonal  matrix  is  stasis  in 
size  class  2  (40-55  mm  clams),  indicating  that  the  difference  in  A  is  mostly  due  to 
differences  in  this  parameter.  This  result  is  in  agreement  with  the  results  of  the 
elasticity  analysis,  although  they  do  not  necessarily  always  agree.  In  the  seasonal 
matrices,  P2  at  clean  sites  were  less  than  0.3  (except  in  winter  when  P2  was  0.9)  while 
it  was  greater  than  0.7  in  all  seasons  at  contaminated  sites.  The  importance  of  this 
element  is  highest  in  summer  and  fall.  Positive  elements  of  are  ones  which  were 
larger  in  ,  and  negative  ones  were  larger  in  A^^ .  Although  elements  of  are 
larger  for  reproduction  and  growth  through  the  size  classes,  changes  in  these  elements 
do  not  contribute  as  much  to  changes  in  A  as  do  the  elements  for  stasis,  which  are 
higher  in  A^^ . 
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log  (Total  PAH) 


Figure  4.8:  Population  Growth  Rate  vs.  PAH  Concentration.  Population  growth 
rates  for  each  site  are  plotted  against  the  logarithm  of  total  PAH  concentration  (in 
ng/g  dry  weight)  at  the  site.  Sites  are,  from  left  to  right,  Wellfleet  Harbor,  Barnstable 
Harbor,  Neponset  River,  Saugus  River,  and  Fort  Point  Channel. 
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Table  4.17:  Annual  Matrices.  Measured  recruitment  from  1995  was  used  to  calculate 
Ri,  then  the  four  seasonal  matrices  were  multiplied  together  for  each  site.  Parameters 
in  these  matrices  include  all  possible  ways  clams  could  make  each  transition  during 
one  year.  Elements  below  the  first  subdiagonal  represent  growth  through  more  than 
one  stage,  which  is  not  possible  in  one  season,  but  common  over  one  year. _ 


Barnstable  Harbor 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.001 

0.002 

0.002 

0.007 

0.004 

0.009 

0.710 

0.079 

0 

0.001 

0.002 

0.002 

0.005 

0.010 

0.286 

0.338 

0.026 

0.038 

0.012 

0 

0 

0 

0 

0 

0.008 

0.027 

0.009 

0 

0 

0 

0 

0 

0.002 

0.006 

0.002 

0 

0 

0 

0 

0 

0 

0 

0.013 

0.049 

0.026 

0.090 

Fort  Point  Channel 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.003 

0.005 

0.008 

0.013 

0.028 

0.0430 

1.000 

0 

0 

0 

0 

0 

0 

0 

0 

0.010 

0.255 

0.503 

0.733 

0 

0 

0 

0 

0 

0 

0 

0.001 

0.005 

0 

0 

0 

0 

0 

0 

0.001 

0.008 

0.004 

0 

0 

0 

0 

0 

0.001 

0.004 

0.005 

0.0165 

Neponset  River 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.034 

0.039 

0.088 

0.100 

0.311 

0.333 

0.833 

0.556 

0.123 

0.113 

0 

0 

0 

0 

0.083 

0.292 

0.799 

0.793 

0.699 

0 

0 

0 

0 

0 

0 

0.008 

0.087 

0.545 

0 

0 

0 

0 

0 

0 

0 

0.005 

0.375 

0 

0 

0 

0 

0 

0 

0.002 

0.125 

0.177 

Saugus  River 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.001 

0.002 

0.004 

0.005 

0.012 

0.022 

1.000 

0.833 

0.139 

0.096 

0 

0 

0 

0 

0 

0.167 

0.695 

0.771 

0.829 

0 

0 

0 

0 

0 

0 

0.020 

0.132 

1.000 

0 

0 

0 

0 

0 

0 

0.001 

0.008 

0.622 

0 

0 

0 

0 

0 

<  0.001 

0.001 

0.292 

0.530 

Wellfleet  Harbor 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1.000 

0 

0 

0 

0 

0 

0 

0 

0 

0.088 

0.001 

<  0.001 

0 

0 

0 

0 

0 

0 

0 

0.004 

0.003 

0 

0 

0 

0 

0 

0 

0.002 

0.008 

0.018 

0 

0 

0 

0 

0 

0 

0.001 

0.009 

0.018 

0.014 

147 


Table  4.18:  Population  Growth  Rates.  Recruits  per  adult,  r,  and  population  growth 
rates,  A,  calculated  using  measured  recruitment  rates  for  each  site,  the  minimum 
measured  recruitment  rate,  r  =  0,  and  the  maximum  measured  recruitment  rate, 
r  =  2.3.  _ 


Site 

r 

^min 

^max 

Barnstable  H. 

0.04 

0.09 

0.09 

0.21 

Wellfleet  H. 

0 

0.01 

0.01 

0.07 

Saugus  R. 

0.09 

1.00 

1.00 

1.08 

Neponset  R. 

2.3 

0.84 

0.70 

0.84 

Fort  Point  C. 

0.21 

0.74 

0.73 

0.76 

4.4  Discussion 

Numerous  studies  on  marine  bivalves  have  shown  negative  effects  of  contaminant 
exposure  on  individual  physiological  processes.  Populations  of  Mya  arenaria  living 
at  contaminated  sites  are  predicted  to  be  composed  of  unhealthy  individuals.  In 
fact,  McDowell  and  Shea  (1996)  report  that  physiological  condition  of  M.  arenaria, 
as  measured  by  digestive  gland-gonad  condition  index,  dry  weight  of  the  digestive 
gland-gonad  complex,  and  the  lipid  weight  of  the  digestive  gland-gonad  complex,  was 
significantly  lower  at  the  contaminated  sites  used  in  this  study  than  at  the  clean  sites. 
These  results  are  similar  to  those  reported  by  Capuzzo  et  al.  (1989)  for  transplanted 
mussels  {Mytilus  edulis)  in  New  Bedford  Harbor,  where  reduced  condition  indices 
and  lipid  accumulation  correlated  with  reduced  reproductive  effort.  Clams  at  the 
contaminated  sites  were  not  only  energetically  inferior  to  clams  at  clean  sites,  but 
also  had  significantly  higher  incidence  of  the  disease  hematopoietic  neoplasia,  and  of 
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Table  4.19:  Life  Table  Response  Experiment  Decomposition  Matrices.  Elements  in 
decomposition  matrices  correspond  to  seasonal  matrix  elements.  The  sum  of  the 
elements  of  all  these  matrices  is  approximately  equal  to  the  difference  between  the 
population  growth  rates  of  the  “clean”  and  “contaminated”  sites.  Negative  values  are 
parameters  that  were  larger  in  the  “clean”  matrices;  positive  values  are  parameters 
that  were  larger  in  the  “contaminated”  matrices. 

Spring  0  0  0  <  0.001  -0.001  -0.001  -0.001  -0.002 

0000  0  0  0  0 

0000  0  0  0  0 

0  0  0  0.003  0  0  0  0 

0  0  0  -0.002  0.095  0  0  0 

0  0  0  0  -0.01  0.02  0  0 

0  0  0  0  -0.002  -0.006  0.007  0 

0  0  0  0  0  -0.005  -0.006  0.001 

Summer  0  0  0  <  0.001  0.002  0.003  0.006  0.013 

0000  0  0  0  0 

0000  0  0  0  0 

0  0  0  0.012  0  0  0  0 

0  0  0  0.019  0.176  0  0  0 

0  0  0  0  -0.001  0.047  0  0 

0  0  0  0  -0.001  0.002  0.013  0 

0  0  0  0  0  -0.002  0.003  0.013 

Fall  0  0  0  0  0  0  0  0 

0000  0  0  0  0 

0000  0  0  0  0 

0  0  0  0.006  0  0  0  0 

0  0  0  -0.002  0.159  0  0  0 

0  0  0  0  -0.004  0.006  0  0 

0  0  0  0  0  -0.003  0.014  0 

0  0  0  0  0  0  <  0.001  0.011 

Winter  0000  0  0  0  0 

0000  0  0  0  0 

0000  0  0  0  0 

0  0  0  0.001  0  0  0  0 

0  0  0  -0.001  0.007  0  0  0 

0  0  0  0  0  0.004  0  0 

0  0  0  0  0  -0.002  -0.002  0 

0  0  0  0  0  0  0.001  -0.002 
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cell  proliferation,  an  abnormality  in  gonadal  tissue.  Despite  the  significance  of  these 
physiological  impairments,  we  found  that  population  growth  rates  were  highest  at  the 
contaminated  sites. 

Higher  population  growth  rates  at  contaminated  sites  are  suggested  by  the  higher 
densities  of  clams  present  at  these  sites.  However,  the  reason  that  these  populations 
are  dense  was  unknown,  and  could  not  be  predicted  by  information  on  individual 
physiological  processes  alone.  By  integrating  clam  survival,  growth  and  reproduction 
in  matrix  population  models,  we  found  that  survival  of  clams  plays  the  largest  role 
in  determining  population  growth  rates.  This  was  confirmed  by  elasticity  analysis 
and  LTRE  decomposition.  Clams  at  contaminated  sites  had  higher  survial  rates  than 
clams  at  clean  sites.  Although  shell  growth  rates  were  higher  at  clean  sites,  survival 
plays  a  larger  role  in  population  growth  rate,  leading  to  faster-growing  populations 
at  the  contaminated  sites.  Fecundities  and  recruitment  did  not  differ  between  con¬ 
taminated  and  clean  sites,  but  even  if  they  had,  the  impact  on  population  growth 
rate  would  be  small  because  of  small  elasticities  to  the  reproductive  parameters  in 
the  model.  Although  many  studies  have  shown  laboratory  effects  of  contaminants  on 
reproduction,  larval  survival,  and  energy  budgets  on  marine  bivalves,  the  strength 
of  these  effects  on  M.  arenaria  and  their  importance  in  the  field  were  insufficient  to 
cause  low  population  growth  rates. 

Sub— lethal  changes  in  clam  physiology  cannot  be  used  to  explain  differences  in 
survival  rates.  Factors  directly  causing  death  during  a  three-month  time  period 
must  be  responsible  for  the  differences  in  survival  between  these  sites.  It  is  possible 
that  different  survival  rates  could  be  due  to  different  predation  rates.  Although 
predators  were  not  carefully  monitored  in  this  study,  we  could  count  mortality  due 
to  moon  snails  {Lunatia  heros)  because  moon  snails  leave  a  characteristic  drill  hole 
in  the  shells  of  clams  that  they  have  eaten.  For  example,  25%  of  summer  mortality 
was  due  to  moon  snail  predation  at  Wellfleet,  but  evidence  of  moon  snail  predation 
was  never  observed  at  the  contaminated  sites.  Clams  at  Wellfleet  Harbor  were  also 
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heavily  preyed  upon  by  the  milky  ribbon  worm,  Cerebratulus  lacteus,  during  the 
spring  and  summer  (pers.  obs.  and  pers.  com.,  Paul  Sommerville,  Wellfleet  Shellfish 
Biologist).  This  worm  appears  sporadically  in  high  abundance  on  clam  flats  and  has 
been  know  to  eliminate  clam  populations  with  low  recruitment  rates  (Rowell  and 
Woo,  1990).  Although  C.  lacteus  is  a  characteristic  member  of  the  fauna  found  in 
contaminated  sites  in  the  New  York  Bight  (Chang  et.  al,  1992),  it  was  not  observed 
at  the  contaminated  sites  in  this  study.  It  is  possible  that  predators  are  more  sensitive 
to  contaminants  than  are  clams.  The  hypothesis  that  high  population  growth  rates  of 
clams  at  Neponset  River  and  Saugus  River  are  due  to  lower  abundance  of  predators 
deserves  further  investigation. 

The  comparison  between  the  Xr  and  Xmin  and  Xmax  shows  that  the  value  used 
for  recruitment  in  the  model  can  be  changed  over  three  orders  of  magnitude  and  the 
population  growth  rate  is  changed  very  little.  Errors  in  estimation  of  recruitment  of 
several  orders  of  magnitude  would  therefore  not  change  the  observed  trends  in  popu¬ 
lation  growth  rates.  This  is  backed  up  by  the  elasticity  analysis.  Several  assumptions 
had  to  be  made  to  estimate  recruitment,  but  since  recruitment  has  so  little  impact 
on  population  growth  rates,  the  validity  of  these  assumptions  is  not  very  important 
to  the  results  of  the  model.  The  fact  that  Barnstable  and  Wellfleet  clams  are  repro¬ 
ducing  in  two  seasons  rather  than  one  apparently  had  little  importance,  since  these 
sites  still  have  dramatically  lower  population  growth  rates  than  the  Boston  sites.  The 
difference  in  spawning  times  is  worthy  of  note,  however  because  populations  of  M.  are- 
naria  in  Maine  are  not  even  spawning  as  late  as  the  Boston  area  sites  (pers.  comm., 
B.  Beale).  Spawning  is  triggered  by  water  temperatures,  so  clams  in  Massachusetts 
should  be  spawning  before  clams  in  Maine,  where  waters  are  colder. 

The  results  reported  here  indicate  that  populations  may  be  decreasing  quite 
rapidly  at  most  of  the  sites,  and  only  growing  at  Saugus.  Population  surveys  made 
in  June  1996  (Figure  4.7),  showed  that  the  projections  of  the  1995  model  were  not 
unreasonable.  A  rapid  decline  in  the  Wellfleet  population  was  expected,  and  in  fact. 
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no  clams  were  found  in  the  study  area  in  June  1996.  The  Barnstable  population  was 
also  projected  to  decrease  rapidly.  In  June  1996,  the  total  population  density  had 
increased  to  9.6  ±  3.41  clams  per  core.  This  was  due  to  a  large,  successful  recruitment 
event,  since  70%  of  the  clams  were  obviously  less  than  a  year  old  (mostly  <30  mm). 
Only  about  30%  of  the  clams  in  1996  were  in  size  classes  3,  4,  and  5.  This  roughly 
translates  to  a  population  density  decrease  to  60%  of  that  of  the  previous  year,  if 
recruitment  is  ignored,  which  is  a  large  decrease,  as  predicted  by  the  model. 

The  low  population  growth  rates  reported  here  are  also  not  totally  accurate,  since 
Ripley  and  Caswell  (in  prep.)  have  shown  that  using  highly  variable  stochastic  recruit¬ 
ment  in  a  matrix  model  instead  of  a  single  value  for  recruitment  allows  populations 
to  grow  at  a  faster  rate,  and  is  a  more  realistic  representation  of  population  processes. 
This  technique  was  not  used  here  because  it  requires  data  on  the  range  of  intensity  in 
recruitment  events  over  many  years,  which  were  unavailable  for  these  sites.  However, 
recruitment  was  shown  to  vary  from  year  to  year  in  this  study.  The  comparison  be¬ 
tween  sites  would  probably  have  the  same  result  using  the  stochastic  method  unless 
extent  of  variability  in  recruitment  differs  substantially  between  sites. 

Populations  of  bivalve  molluscs  chronically  exposed  to  contaminated  habitats  may 
be  highly  resistant  to  contaminant  effects  (McDowell  Capuzzo,  1996).  In  this  study, 
we  found  that  population  growth  of  Mya  arenaria  was  enhanced  by  contamination. 
It  is  possible  that  genetic  differences  between  populations  at  these  sites  contribute 
to  the  observed  contaminant  tolerances.  Circulation  studies  (Geyer  et  ai,  1992)  of 
the  Massachusetts  Bay  system  suggest  that  retention  time  in  the  Bay  is  on  the  or¬ 
der  of  three  weeks — longer  than  the  two-three  week  larval  period  of  M.  arenaria. 
This  suggests  that  larvae  released  in  Cape  Cod  Bay  are  probably  retained  there  until 
they  settle.  Boston  area  populations  are  probably  subsidized  by  larvae  from  popu¬ 
lations  farther  north.  Their  larvae  are  probably  swept  southward  by  the  consistent 
counterclockwise  circulation  pattern  in  Massachusetts  Bay  (Geyer  et  al,  1992).  All 
populations  probably  do  not  come  from  one  well-mixed  larval  pool,  although  Boston 
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sites  and  Cape  Cod  sites  may  see  separate,  mixed  pools.  These  conclusions  are  based 
on  the  assumption  that  larvae  released  over  tidal  mudflats  become  entrained  into  the 
general  circulation  pattern  of  deeper  water. 

The  supply  of  larvae  imported  to  contaminated  sites  from  sites  farther  north 
may  be  sufficient  to  overcome  any  potential  reproductive  impairments.  In  a  larval 
settlement  assay,  Wintermeyer  et  al.  (1996)  showed  that  M.  arenaria  settled  equally 
well  on  sediment  from  all  of  the  sites  studied  here,  except  for  Wellfleet  Harbor.  High 
organic  matter  (16.8%)  and  low  oxygen  levels  in  the  Wellfleet  sediments  may  explain 
the  reduced  settlement  success  at  this  site.  This  suggests  that  although  larvae  have 
been  shown  to  be  highly  susceptible  to  contaminants  in  laboratory  experiments,  they 
are  able  to  settle  and  metamorphose  when  transported  to  contaminated  sites  under 
field  conditions. 

Although  many  studies  have  shown  deleterious  effects  on  contaminant  exposure 
on  marine  bivalves,  the  kind  of  sublethal  impairments  demonstrated  by  these  studies 
were  not  observed  to  lead  to  reduced  population  levels  at  contaminated  sites  in  this 
study.  My  a  arenaria  is  a  species  which  is  quite  tolearant  to  contaminant  stress,  so  our 
results  may  not  apply  to  less  hardy  species,  however  it  is  impossible  to  study  natural 
populations  of  those  species  at  contaminated  sites  in  the  field.  We  have  shown  that 
it  is  important  to  consider  other  ecological  factors,  such  as  predation  rates,  and  life 
cycle  factors,  such  as  potential  variation  in  recruitment,  when  evaluating  contaminant 
effects  on  field  populations. 

4.5  Summary 

1.  Clams  at  clean  sites  grew  (in  shell  length)  significantly  faster  than  those  at  contam¬ 
inated  sites,  but  clams  at  contaminated  sites  had  higher  survival  rates.  No  differences 
were  found  in  fecundity  or  recruitment  (by  area),  although  reproductive  seeison  was 
shorter  at  contaminated  sites  and  recruitment  (per  adult)  did  differ. 
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2.  Population  growth  rates  calculated  in  matrix  population  models  were  highest  at 
contaminated  sites. 

3.  Differences  between  mortality  rates  (possibly  due  to  predation)  and  long-term 
recruitment  patterns  at  the  sites  may  have  more  of  an  influence  on  the  dynamics  of 
clam  populations  than  contaminant  exposure. 
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Chapter  5 


Conclusion 


The  results  of  this  thesis  show  that  the  life  history  traits  of  a  marine  bivalve  can 
have  important  consequences  for  the  effects  changes  in  the  environment  will  have  on 
populations  of  the  bivalve.  For  the  soft  shell  clam,  a  broadcast  spawner,  the  highest 
elasticities  are  for  stasis,  that  is,  survival  without  growth.  According  to  the  literature 
data  on  several  broadcast  spawning  species,  this  was  especially  true  for  those  species 
with  longer  life  spans.  Life  span  is  important  because  a  broadcast  spawning  reproduc¬ 
tive  mode  leads  to  uncertain  recruitment  success  in  each  year.  However,  the  longer 
the  life  span,  the  more  chances  the  population  has  to  experience  a  successful  recruit¬ 
ment  event.  In  fact,  variation  in  recruitment  levels  actually  increases  the  population 
growth  rate  over  that  at  mean  recruitment  levels.  Analysis  of  data  on  vital  rates  of 
clams  at  clean  and  contaminated  sites  showed  that  in  fact,  survival  rates  were  the 
most  important  contributors  to  differences  between  the  sites  studied.  Although  clams 
at  clean  sites  grew  (in  shell  length)  faster,  the  populations  at  contaminated  sites  were 
growing  faster  because  of  significantly  higher  survival  rates.  Lower  predation  rates  at 
contaminated  sites  are  a  possible  cause  of  the  different  survival  rates.  Populations  of 
clams  at  contaminated  sites  appear  to  be  structured  by  the  combination  of  long-term 
recruitment  patterns  and  other  ecological  factors  such  as  predation,  rather  than  only 
by  contaminant  effects  on  individual  clams.  A  different  study  animal  with  different 
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life  history  traits,  such  as  a  brooding  bivalve,  might  show  a  more  direct  response  to 
contaminant  exposure  if  the  effects  of  the  toxicant  were  on  a  process  that  contributed 
strongly  to  population  growth  rate. 
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Appendix  A 


Untangling  Microhabitat 
Variability  in  Growth  Rates  of  Soft 
Shell  Clams  from  Contamination 
Affects  on  Growth 

A.l  Introduction 

Determining  the  factors  responsible  for  patterns  in  growth  rates  of  bivalves  over  large 
and  small  spatial  scales  has  been  difficult  because  variability  is  high  and  multiple 
factors  contribute  to  growth  rates.  Brousseau  failed  to  find  any  latitudinal  gradient 
in  growth  rates  of  the  ribbed  mussel,  Geukensia  demissa  (Brousseau,  1984)  or  the 
soft-shell  clam,  Mya  arenaria  (Brousseau,  1979),  concluding  that  local  environmental 
conditions  play  a  large  role  in  bivalve  growth.  However,  in  a  more  ambitious  study  of 
M.  arenaria  growth  from  Maryland  to  Nova  Scotia,  a  trend  of  faster  growth  at  lower 
latitudes  was  detected  using  principle  component  analysis  (Appledoorn,  1983).  The 
two  other  components  explaining  variation  were  sediment  siltiness  and  hydrocarbons, 
both  negatively  correlated  to  growth  rates.  Growth  rate  depressions  in  G.  demissa  in 
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central  Jamaica  Bay,  New  York  (near  Kennedy  airport)  have  been  suggested  to  be  in 
part  due  contaminant  stress  (Franz  and  Tanacredi,  1993).  Hydrocarbon  contamina¬ 
tion  has  also  been  shown  to  reduce  growth  rates  of  clams  transplanted  to  oiled  sites, 
(Dow,  1975;  Appledoorn,  1981)  suggesting  that  impaired  growth  may  be  an  indicator 
of  chronic  contaminant  exposure. 

Untangling  any  contaminant  effect  on  growth  rate  from  a  local  site  effect  requires 
that  other  factors  known  to  modify  growth  rates  also  be  controlled,  or  at  least  mea¬ 
sured.  Early  researchers  suggested  that  environmental  factors  such  as  food  supply 
(Kellogg,  1905),  current  speed  (Belding,  1930),  and  sediment  qualities  (Newcombe, 
1935)  were  responsible  for  differences  in  growth  rates  between  sites.  Spear  and  Glude 
(1956)  demonstrated  in  a  reciprocal  transplant  experiment  that  environmental  qual¬ 
ities  of  sites,  rather  than  genetic  differences  between  M.  arenaria  sub-populations, 
were  responsible  for  differing  growth  rates.  That  environment,  not  heredity,  controls 
growth  rates  has  been  supported  more  recently  for  Mytilus  edulis  (Kautsky  et  al, 
1990),  and  for  Mercenaria  mercenaria  (Landry  et  al,  1993).  A  more  sophisticated 
study  on  M.  mercenaria  (Rawson  and  Hilbish,  1991)  detected  heritable  variation  in 
juvenile  growth  rates,  but  only  at  the  habitats  where  high  and  low  extremes  in  growth 
were  measured.  Overall,  environmental  differences  accounted  for  more  than  50%  of 
the  variation  in  growth  rates  in  this  study. 

As  an  infaunal,  intertidal  filter  feeder,  environmental  factors  likely  to  affect  Mya 
arenaria  growth  rates  include  tidal  elevation,  phytoplankton  abundance,  flow  rates, 
and  sediment  composition.  Tidal  elevation  is  important  because  animals  living  at 
higher  tidal  levels  will  have  a  shorter  period  of  time  each  day  in  which  they  can  Alter 
food,  respire,  and  flush  wastes.  Newcombe  (1935)  found  that  M.  arenaria  in  the  Bay 
of  Fundy  at  the  eight  foot  tidal  level  grew  50%  more  per  year  than  clams  at  the 
sixteen  foot  tidal  level,  and  that  variation  in  growth  was  directly  correlated  with  time 
of  submergence.  Tidal  height  is  also  known  to  influence  growth  rates  in  G.  demissa 
(Stiven  and  Gardner,  1992),  and  most  of  the  difference  is  explained  by  difference 
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in  submergence  time  (Lin,  1989).  A  mathematical  model  for  M.  edulis  populations 
that  incorporated  food  density  and  time  of  exposure  to  air  showed  that  differences 
in  growth  rate  were  due  to  food  availability,  not  physiological  impairments  due  to 
exposure  (Ross  and  Nisbet,  1990). 

Current  flow  was  cited  by  Belding  (1930)  as  the  most  important  factor  in  M. 
arenaria  growth.  More  recent  work  has  explicated  this  phenomenon.  For  example, 
where  higher  flow  rates  deliver  a  higher  total  amount  of  food  to  the  clams,  faster  flows 
lead  to  faster  growth  rates  {Mya:  Emerson,  1990;  Mercenaria:  Grizzle  and  Morin, 
1989).  Judge  et  al.  (1992)  demonstrated  by  altering  rates  of  flow  at  contiguous 
locations  at  one  site  that  current  speed  alone  does  not  lead  to  different  growth  rates 
of  M.  mercenaria  if  chlorophyll  a  concentrations  are  the  same.  Thus  both  current 
flow  and  tidal  height  are  important  in  the  amount  of  food  available  to  the  animal  at 
the  location,  not  in-and-of  themselves. 

Sediment  composition  has  also  been  shown  to  affect  growth  rates  in  M.  arenaria. 
Newcombe  (1935)  showed  that  clams  in  muddy  soils  with  a  silty  surface  layer  grew 
more  slowly  than  in  sandier  soils.  In  a  sediment  and  clam  reciprocal  transplant 
experiment.  Swan  (1952)  observed  that  sediment  grain  size  may  alter  growth  rate. 
In  coarser  sediments,  clams  had  to  grow  thicker  shells  to  overcome  stronger  shell 
abrasion,  and  could  not  grow  as  fast  as  a  result.  Unfortunately  this  study  had  very 
small  sample  sizes  and  the  results  cannot  be  considered  conclusive.  Sediment  type  is 
a  function  of  flow  rates,  which  are  also  known  to  influence  growth  rates,  so  decoupling 
the  effect  of  sediment  type  from  that  of  flow  rate  is  challenging. 

Resuspention  of  sediment  has  been  reported  to  have  both  a  positive  and  a  neg¬ 
ative  effect  on  growth  rates.  Under  turbid  conditions,  clams  may  reduce  filtration 
rate  and  therefore  food  intake,  reducing  growth  potential  (Grant  and  Thorpe,  1991). 
Disturbances  that  resuspend  benthic  detritus  or  diatoms  may  provide  additional  food 
to  clams,  promoting  growth  (Emerson,  1990).  The  relative  importance  of  these  effects 
is  yet  to  be  determined. 
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Food  availability  has  been  found  to  be  the  strongest  correlate  to  growth  rates 
in  many  bivalves.  A  series  of  papers  (summarized  in  Thompson  and  Nichols,  1988) 
on  Macoma  balthica,  which  co-occurs  with  M.  arenaria,  have  shown  that  although 
temperature  is  strongly  correlated  with  growth  rates,  it  is  really  phytoplankton  pro¬ 
duction  (which  is  related  to  temperature)  that  controls  growth  rates.  Chlorophyll 
a  concentrations  in  the  water  above  clam  flats  have  been  found  to  correlate  most 
strongly  to  differences  in  growth  rates  in  M.  mercenaria  (Landry  et  al,  1993)  and 
Yoldia  notabilis,  a  subtidal,  infaunal,  deposit  feeding  clam  (Nakaoka,  1992). 

Several  environmental  factors  have  been  shown  to  not  be  important  in  determining 
growth  rates  (when  making  comparisons  within  seasons).  Newcombe  (1935)  ruled  out 
temperature,  salinity  and  sulfide  concentrations.  Temperature  and  salinity  were  also 
found  to  not  have  an  effect  on  growth  rates  in  Yoldia  notabilis  (Nakaoka,  1992). 
Although  clam  density  can  affect  growth  rates,  particularly  of  younger  bivalves  {e.g. 
Goshima,  1982;  Stiven  and  Gardner,  1992),  densities  can  be  controlled  in  experimental 
design.  It  has  also  been  shown  in  many  studies  that  growth  rates  decrease  with 
increasing  size  {e.g.  Ripley  and  McDowell,  in  prep.),  but  using  clams  in  a  limited  size 
range  controls  this  factor. 

A  synthesis  of  the  importance  of  these  various  factors  was  reported  by  Grizzle  and 
Morin  (1989).  In  a  factorial  design  field  experiment,  they  tested  the  interaction  of 
tidal  currents,  seston,  and  sediment  type.  Horizontal  seston  fluxes  were  found  to  be 
the  significant  factor.  In  summary,  the  potential  environmental  difference  which  must 
be  controlled  for  in  order  to  evaluate  growth  rates  is  the  amount  of  food  available. 

A  recent  study  of  populations  of  M.  arenaria,  around  Massachusetts  Bay  found 
significantly  higher  growth  rates  at  uncontaminated  sites.  (Ripley  and  McDowell,  in 
prep.).  However,  in  that  study  growth  measurements  were  made  at  only  one  small 
region  of  the  mud  flat  at  each  site,  and  other  factors  known  to  influence  growth  were 
not  measured.  Environmental  factors  unrelated  to  contaminant  concentrations  may 
have  caused  the  observed  growth  rates.  The  study  presented  here  was  done  to  de- 
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termine  whether  within-site  variation  in  growth  rates  was  as  large  as  the  observed 
between-site  difference  in  growth  rates,  when  other  environmental  factors  were  more 
carefully  controlled.  Growth  rates  over  three  months  in  the  summer  were  measured 
for  small  M.  arenaria  at  five  sub-sites  in  one  contaminated  and  one  uncontaminated 
site.  Sediments  were  analyzed  for  total  organic  carbon  as  a  measure  of  hydrocarbon 
load,  and  nitrogen  as  a  measure  of  phytoplankton  detritus  (food  supply).  Mud  frac¬ 
tion  of  sediment  was  measured,  as  was  immersion  time,  since  these  factors  may  also 
contribute  to  growth  rates. 

A. 2  Methods 

A. 2.1  Sites 

Two  sites  were  chosen  for  this  study  based  on  information  on  sediment  contaminant 
levels  (McDowell  and  Shea,  1997).  The  contaminated  site,  the  Neponset  River  in 
Quincy,  MA,  has  a  total  sediment  PAH  concentration  of  about  1500  ng/g  dry  weight, 
while  the  uncontaminated  site,  Barnstable  Harbor,  MA,  has  a  total  sediment  PAH 
concentration  of  about  350  ng/g  dry  weight.  Five  locations  were  selected  randomly 
(with  respect  to  sediment  qualities)  at  each  site  along  the  shoreline  at  approximately 
equal  tidal  heights  (exposed  at  mean  low  water)  as  sub-sites  or  plots. 

A. 2. 2  Clam  Growth 

The  smallest  clams  present  in  the  populations  at  the  two  sites  were  collected  in 
June  1997.  Two  hundred  forty  clams  from  each  site  were  brought  to  the  laboratory 
overnight  and  individually  numbered  in  indelible  ink  and  measured.  For  each  site, 
clams  were  randomly  divided  into  groups  of  twelve.  Each  group  was  planted  in  one 
of  four  sediment-filled  mesh  bags  at  each  plot  (see  Ripley  and  McDowell,  in  prep, 
for  more  detailed  methods).  Bags  were  collected  after  three  months,  and  clams  were 
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sorted  from  the  sediment  by  hand  and  measured. 


A. 2. 3  Immersion  Time 

Immersion  time  was  measured  over  one  lunar  month  at  each  plot.  Measurements  were 
made  using  simple  instruments  constructed  by  B.  Lancaster  at  Woods  Hole  Oceano¬ 
graphic  Institution,  using  Omron  electronic  hour  meters.  Meters  were  sealed  into 
PVC  pipes  such  that  coverage  by  seawater  would  close  the  circuit  between  protrud¬ 
ing  electrodes  and  start  the  meter.  Instruments  were  tested  in  the  laboratory  before 
deployment,  and  they  all  counted  exactly  the  same  times.  Instruments  were  deployed 
by  burying  them  into  the  sediment  with  just  the  electrodes  sticking  up  above  the 
sediment  surface,  and  attaching  each  one  to  a  stake  with  tie-wraps.  Previous  mea¬ 
surements  of  tidal  immersion  between  sites  showed  less  than  10%  variation  between 
monthly  readings,  so  only  one  measurement  was  made. 

A. 2. 4  Semiment  Samples 

When  plots  were  set  up  in  June,  1997,  duplicate  sediment  samples  were  taken  at  each 
plot.  Sediment  from  the  surface  5  cm  was  scooped  into  solvent-rinsed  glass  jars  with 
a  solvent-rinsed  aluminum  spatula.  Sediment  was  frozen  at  -10°C  immediatley  upon 
return  to  the  lab.  Prior  to  analyses,  sediments  were  refrigerated  at  4°C  to  thaw  on 
November  5,  1997.  After  thawing,  each  sample  was  homogenized  by  stirring  with  a 
solvent-rinsed  aluminum  spatula.  Aliquots  of  the  sample  were  taken  for  the  following 
analyses. 

Total  Organic  Carbon  and  Nitrogen 

Approximately  2  g  of  sediment  was  dried  for  24  hours  at  60°  C.  Dry  sediment  was 
ground  with  a  mortar  and  pestle  and  approximately  1  g  was  placed  in  a  plastic  urine 
sample  cup.  Inorganic  carbon  was  dissolved  by  adding  20  ml  1  N  acetic  acid  for  24 
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hours.  Up  to  another  20  ml  of  acetic  acid  were  added  until  no  more  bubbles  were 
observed.  Sediments  were  rinsed  three  times  with  distilled  water,  and  filtered  onto 
0.45  jira.  millipore  paper.  Samples  were  dried  for  24  hours  at  60°C,  weighed,  and 
ground  with  a  mortar  and  pestle.  C:H:N  analysis  was  performed  on  an  Eager  200 
automated  analyzer. 

Sediment  Grain  Size 

Approximately  20  g  sediment  were  placed  in  tared  plastic  weigh  boats,  weighed, 
and  dried  for  48  hours  at  60°C.  The  dry  sediment  was  weighed,  then  rinsed  with 
approximately  250  ml  water  in  a  63/rm  seive.  Sediment  remaining  in  the  seive  was 
rinsed  onto  a  pre-weighed  glass  fiber  filter.  Sediment  on  filter  (and  filter)  were  replaced 
into  the  weigh  boat  and  dried  for  48  hours  at  60°C.  Fraction  of  sample  less  than  63/im 
(the  mud  fraction)  was  calculated  as  1  —  (remaining  sediment) /(start  weight). 

A.  3  Results 

A. 3.1  Clam  Growth 

In  June,  the  mean  length  of  clams  deployed  in  Barnstable  Harbor  was  27.9  ±  0.23 
(one  standard  error),  and  the  mean  length  of  clams  in  Neponset  River  was  37.5  ± 
0.17.  In  September,  46  clams  were  found  alive  from  Barnstable,  and  91  were  found 
alive  in  Neponset.  At  least  one  clam  was  found  at  each  plot  at  both  sites.  At 
the  start  of  the  study  clam  sizes  were  different  enough  between  sites  that  growth 
rates  might  be  different  simply  due  to  starting  sizes.  Clam  growth  rates  within  each 
site  between  plots  were  compared  by  one-way  ANOVA.  Although  growth  differed 
significantly  between  plots  at  each  site  (Table  A.l),  there  appears  to  be  a  difference 
between  growth  rates  between  sites  as  well  (Figure  A.l).  However,  the  different 
starting  sizes  of  clams  makes  this  apparent  difference  difficult  to  interpret  and  test. 
The  slopes  of  the  regressions  of  growth  on  start  length  (Figure  A. 2)  are  significantly 
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Table  A.l:  ANOVA  table  for  growth  between  plots  for  each  site.  Source  of  variation 
(Source)  is  explained  by  the  sum  of  squares  (SS),  degrees  of  freedom  (df),  mean 
square  (MS),  F  statistic  (F),  and  significance  (p). 


Source 

SS 

df 

MS  F 

P 

Barnstable  Harbor 

Between  plots 

341.5 

4 

85.4  3.08 

0.03 

Residual 

1165 

42 

27.7 

Total 

1506.8 

46 

32.8 

Neponset  River 

Between  plots 

244.5 

4 

61.1  9.9 

<0.001 

Residual 

540.1 

88 

6.1 

Total 

784.6 

92 

8.5 

169 


BH  QN 

Location 

Figure  A.l:  Clam  growth  by  plot.  Mean  growth  rates  (mm/3  mo)  over  all  clams 
recovered  alive  in  September  from  each  plot  (1-5,  from  left  to  right)  for  Barnstable 
Harbor  (BH)  and  Neponset  River  (QN).  Error  bars  represent  one  standard  error  and 
where  they  are  not  shown  the  bar  represents  only  one  clam. 
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Table  A.  2:  ANOVA  table  for  slope  of  regressions  of  growth  on  start  length  between 
sites.  Source  of  variation  (Source)  is  explained  by  the  sum  of  squares  (SS),  degrees 
of  freedom  (df),  mean  square  (MS),  F  statistic  (F),  and  significance  (p). 


Source 

SS 

df 

MS 

F 

P 

Within-|-residual 

1666.9 

136 

12.3 

Site 

365.9 

1 

365.9 

29.8 

<0.001 

Start  length 

375.7 

1 

375.7 

30.6 

<0.001 

Site  X  start  length 

189.7 

1 

189.7 

15.5 

<0.001 

different  (Table  A. 2)  demonstrating  that  clams  differ  in  their  size-dependent  growth 
rates  depending  on  their  location,  so  an  anlysis  of  covariance  cannot  be  done. 

A. 3. 2  Immersion  Time 

Measurements  of  hours  wet  per  day  were  made  once  for  each  plot  at  Neponset  River. 
All  but  one  of  the  tidal  meters  deployed  in  Barnstable  Harbor  were  either  lost  or 
leaked  and  malfuctioned.  Data  collected  previously  (1-5/97)  at  one  of  the  plots  at 
Barnstable  can  be  used,  so  we  have  data  for  two  plots  at  Barnstable  Harbor,  specif¬ 
ically,  20.7  (plot  1)  and  11.2  (plot  4)  hours  wet/day  (mean=15.95,  variance=45.12). 
Water  coverages  at  Neponset  were  21.7,  18.9,  22.9,  17.2,  and  18.3  (plots  1-5)  hours 
per  day  (mean=19.8,  variance=5.76).  A  t-test  comparing  means  at  the  two  sites  was 
not  significant  (p=0.27). 
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Growth 


Start  Length 


Figure  A.2:  Growth  by  start  length.  Growth  (mm/3  mo)  for  each  clam  is  plotted 
against  its  length  at  the  start  of  the  study.  Triangles  are  clams  from  Barnstable 
Harbor  and  squares  are  clams  from  Neponset  River. 
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Table  A. 3:  ANOVA  table  for  mud  fraction  of  sediments.  Source  of  variation  (Source) 
is  explained  by  the  sum  of  squares  (SS),  degrees  of  freedom  (d/),  mean  square  (MS), 
F  statistic  (F),  and  significance  (p). 


Source 

SS 

df 

MS 

F 

P 

Among  sites 

0.001 

1 

0.001 

0.111 

>0.05 

Among  plots  within  sites 

0.072 

8 

0.009 

45 

<0.001 

Within  plots 

0.002 

10 

0.0002 

Total 

0.075 

19 

A. 3, 3  Sediment  Samples 

Sediment  Grain  Size 

Mud  fractions  of  sediments  varied  between  plots,  but  about  the  same  range  (6%  to 
23%)  was  represented  at  each  site  (Figure  A.3).  Mud  fractions  were  compared  with  a 
nested  design  ANOVA  (Table  A.3).  No  significant  variation  was  found  between  sites 
(p>0.05),  but  there  was  a  highly  significant  (p<0.001)  added  variance  component 
among  plots. 

Sediment  Organic  Carbon  and  Nitrogen 

Sediments  were  about  1%  organic  carbon  (by  weight)  and  about  0.1%  nitrogen.  The 
two  plots  (4  and  6)  with  high  levels  of  carbon  and  nitrogen  (Figure  A. 4)  are  also  high 
in  mud  content  (Figure  A.3).  Although  there  was  slightly  more  carbon  and  nitrogen 
at  Neponset  than  at  Barnstable,  the  difference  between  sites  was  not  significant  for 
either  factor  when  tested  in  a  nested  design  ANOVA  (Table  A. 4).  There  was  a 


173 


Location 


Figure  A.3:  Sediment  Grain  Size.  Bars  show  the  mean  of  two  replicate  samples  of 
mud  fraction  (g  <63/xm  per  gram  sediment)  for  the  five  plots  (1-5,  left  to  right)  at 
Barnstable  (BH)  and  Neponset  (QN).  Error  bars  are  one  standard  error. 


Plot 


Figure  A. 4:  Sediment  organic  carbon  and  nitrogen.  Organic  carbon  (g/g  sediment) 
is  shown  in  the  top  panel,  and  nitrogen  (g/g  sediment)  is  shown  in  the  lower  panel. 
Points  shown  are  each  of  three  replicate  C:H:N  anlyses  on  each  sediment  sample. 
Filled  squares  are  one  replicate  sediment  sample  from  a  plot,  and  open  circles  in  the 
same  column  are  the  other  replicate  from  the  same  plot.  Plots  are  numbered  1-10, 
where  1-5  are  at  Barnstable  Harbor  and  6-10  are  at  Neponset  River. 
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significant  added  variance  component  for  carbon  due  to  both  plots  (p<0.01)  and 
replicates  (p<0.001),  but  for  nitrogen  these  factors  were  not  significant  (p>0.05). 

A. 4  Discussion 

The  question  this  study  was  designed  to  answer  was  whether  observed  differences  in 
clam  growth  rates  between  contaminated  and  uncontaminated  sites  were  likely  due 
to  contaminant  levels,  or  whether  they  could  be  explained  by  some  of  the  many  other 
physical  factors  that  could  differ  between  sites.  The  results  of  this  study  reinforced 
our  previous  results  (Ripley  and  McDowell,  in  prep.)  in  that  clam  growth  rates  at 
Barnstable  Harbor,  the  uncontaminated  site,  were  higher  than  at  Neponset  River, 
the  contaminated  site.  The  significance  of  this  difference  could  not  be  evaluated  due 
to  differences  in  size  of  clams  used  at  each  site,  however,  in  the  size  range  that  was 
represented  at  both  sites  (30-35  mm),  clams  at  Barnstable  Harbor  grew  about  10  mm 
more  than  clams  at  Neponset  River  did  in  the  three  month  period.  This  difference  in 
growth  could  mean  that  clams  at  Barnstable  reached  reproductive  size  (ca.  40  mm: 
Belding,  1930)  a  year  before  Neponset  River  clams  and  is  probably  an  important 
difference,  regardless  of  statistical  significance. 

Evaluating  whether  the  observed  difference  in  growth  rates  between  sites  could 
be  due  to  any  factor  known  to  influence  growth  rates  in  bivalves  besides  contaminant 
levels  was  the  real  object  of  this  study.  These  factors  are  sediment  composition 
(Appledoorn,  1981),  submergence  time  (Lin,  1989),  and  food  availability  (Judge  et 
al.,  1992).  We  showed  that  mud  fractions  of  sediments  varied  among  plots  at  both 
sites,  but  did  not  differ  between  sites.  Sediment  composition  thus  does  not  explain 
the  observed  pattern  of  clam  growth  rates.  Submergence  times  did  not  differ  between 
sites,  either.  Although  data  were  incomplete  at  Barnstable  Harbor,  missing  data 
points  would  all  have  to  be  close  to  the  lower  value  to  change  the  significance  of 
the  test.  This  would  mean  that  submergence  times  would  be  significanlty  lower  at 
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Table  A. 4:  ANOVA  table  for  organic  carbon  and  nitrogen  in  sediments.  Source  of 
variation  (Source)  is  explained  by  the  sum  of  squares  (SS),  degrees  of  freedom  (df), 
mean  square  (MS),  F  statistic  (F),  and  significance  (p). 


Source 

SS 

df  MS 

F 

P 

Organic  Carbon 

Among  sites 

1 

1.44x10-®  1.44x10-® 

0.36 

>0.05 

Among  plots  w/in  sites 

8 

3.08x10-^  3.85x10-® 

5.0 

<0.01 

Among  reps  w/in  plots 

10 

7.7x10-®  7.7x10-® 

7.0 

<0.001 

W/in  reps 

40 

4.4x10-®  1.1x10-® 

Total 

59 

4.43x10-'* 

Nitrogen 

Among  sites 

1 

1.06x10-® 

1.06x10-® 

2.97 

>0.05 

Among  plots  w/in  sites 

8 

2.86x10-® 

3.57x10-^ 

2.13 

>0.05 

Among  reps  w/in  plots 

10 

1.68x10-® 

1.68x10-^ 

67.2 

>0.05 

W/in  reps 

40 

1.0x10-^ 

2.5x10-^ 

Total 

59 

5.69x10-® 
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Barnstable  Harbor,  which  does  not  explain  the  pattern  of  growth  rates,  because  lower 
submergence  times  should  lead  to  lower  growth  rates  (Lin,  1989). 

Sediment  nitrogen  levels  were  measured  as  an  approximation  to  food  content  of  the 
overlying  water.  No  differences  were  found  between  sites,  or  between  plots,  demon¬ 
strating  that  differences  in  growth  would  not  be  correlated  to  differences  in  sediment 
nitrogen.  Sediment  nitrogen  may  not  be  a  good  enough  aproximation  to  food  sup¬ 
ply  to  the  clams.  Measuring  actual  chlorophyll  a  content  in  water  overlying  mud 
flats  has  been  shown  to  explain  growth  patterns  (Grizzle  and  Morin,  1989),  but  this 
method  was  beyond  the  scope  of  our  interest  in  this  question.  Similarly,  total  organic 
carbon  content  of  sediments  is  only  an  approximate  measure  of  lipophilic  organic 
contaminant  concentrations.  Carbon  was  not  shown  to  differ  between  sites,  but  high 
variation  between  sub-samples  could  have  made  site  differences  undetectable.  How¬ 
ever,  according  to  our  measurements,  organic  carbon  could  not  explain  differences  in 
growth  rates,  either. 

A  factor  we  did  not  evaluate  that  could  play  a  role  in  growth  rates  are  genetic 
differences  between  sites  (Rawson  and  Hilbish,  1991).  It  is  possible  that  clams  at 
Barnstable  are  a  faster-growing  strain,  and  it  would  be  interesting  to  investigate 
this  question  further.  It  is  also  possible  that  some  other  environmental  factor  that 
we  did  not  measure  is  responsible  for  determining  clam  growth  rates,  although  all 
factors  previously  demonstrated  to  be  important  to  bivalve  growth  were  measured  in 
some  way.  The  methods  we  used  to  measure  clam  food  supply  and  sediment  organic 
contaminants  were  only  approximations,  however,  and  more  conclusive  results  might 
be  found  using  more  accurate  techniques. 

A. 5  Summary 

None  of  the  measured  environmental  factors  (sediment  organic  carbon,  sediment  ni¬ 
trogen,  sediment  mud  fraction,  immersion  time)  showed  the  same  pattern  of  consis- 
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tently  higher  values  at  Barnstable  Harbor  plots  as  did  clam  growth  rates. 
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